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Abstract 

When  conducting  a  space  proximity  operation,  developing  high-fidelity  esti¬ 
mates  of  the  dynamical  and  physical  properties  of  a  Resident  Space  Object  (RSO) 
based  on  post-rendezvous  observational  data  acquired,  is  imperative  for  the  under¬ 
standing  of  the  RSO  itself  and  the  operating  environment.  This  research  investigates 
the  estimation  of  relative  motion  dynamics,  rotational  dynamics,  and  the  feasibility 
of  estimating  the  moments  of  inertia  of  a  RSO.  Using  the  Hill-Clohessy-Wiltshire 
equations,  rigid-body  dynamics,  and  estimation  theory,  a  nonlinear  least  squares 
estimation  algorithm  is  implemented  in  the  processing  of  range  data  from  tracked 
observation  points  on  the  RSO  body.  Through  simulation,  it  was  determined  that  ac¬ 
curately  estimating  the  relative  motion  and  rotational  dynamics  is  possible.  However 
directly  estimating  the  moments  of  inertia  using  range  data  proved  to  be  problematic 
and  exposed  a  possible  observability  limitation.  Yet  in  general,  the  solutions  were 
heavily  dependent  on  the  quality  of  the  a  priori  knowledge  as  well  as  the  reduction 
of  solution  ambiguity  through  the  use  of  multiple  observational  data  sets. 
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SPACECRAFT  PROXIMITY  OPERATIONS  USED  TO 


ESTIMATE  THE  DYNAMICAL  &  PHYSICAL  PROPERTIES  OF 
A  RESIDENT  SPACE  OBJECT 


I.  Introduction 

Major  advances  in  the  area  of  relative  navigation  and  rendezvous  are  leading  to  new 
applications  of  spacecraft  missions  which  will  greatly  further  our  capability  to  operate 
in  space.  With  the  developing  engineering  and  technology  of  autonomous  rendezvous 
missions,  the  next  logical  step  is  conducting  autonomous  proximity  operations  on  or 
about  a  Resident  Space  Object  (RSO).  This  is  a  very  exciting  and  promising  field 
with  incredible  futuristic  applications  for  both  the  civilian  and  the  military  space 
communities.  Conducting  a  proximity  operation  requires  the  performing  spacecraft 
to  carefully  function  in  a  very  active  and  dynamic  environment;  where  typically, 
very  little  is  known  for  certain.  The  complexity  of  the  operating  environment  can 
be  successfully  managed  if  the  dynamical  parameters  required  to  operate  near  the 
RSO  can  be  accurately  characterized  and  understood  [1].  Whether  the  proximity 
mission  is  solely  limited  to  making  observations  of  the  RSO  or  if  the  observation 
activity  is  a  precursor  to  a  more  involved  and  intrusive  stage  of  a  proximity  mission, 
understanding  the  dynamics  and  physical  properties  of  the  RSO  provides  valuable 
insight  into  the  operating  environment,  as  well  as  the  RSO  itself. 

This  thesis  will  explore  the  use  of  nonlinear  least  squares  theory  to  estimate 
the  dynamical  and  physical  properties  of  a  RSO.  Specifically,  the  relative  motion  dy¬ 
namics,  rotational  dynamics,  and  moments  of  inertia  will  be  estimated  using  range 
observations  made  by  an  observing  spacecraft  performing  a  space  proximity  opera¬ 
tion  in  the  vicinity  of  the  RSO. 


1-1 


1.1  A  Historical  Perspective 

Several  milestone  space  programs  in  the  history  of  astronautics  have  led  to 
today’s  advances  in  autonomous  rendezvous  and  ultimately  to  the  feasibility  of  per¬ 
forming  autonomous  proximity  operations  and  observations.  First  demonstrating  the 
initial  concepts  of  rendezvous  and  docking  techniques  were  the  Gemini  and  Apollo 
programs  [2],  The  Space  Shuttle  program  then  furthered  the  versatility  and  fre¬ 
quency  of  rendezvous  and  proximity  operation  missions  in  low  earth  orbit  (LEO) 
[2],  The  move  from  a  human-controlled  operation  to  a  new  era  of  autonomous  ren¬ 
dezvous  and  autonomous  proximity  operations  has  been  underway  as  seen  in  recent 
programs  such  as  the  XSS-11  mission  and  the  Hayabusa  mission  [2,  3]. 

1.1.1  The  Gemini  &  Apollo  Programs.  In  1965,  a  modified  Titan  II  In¬ 
tercontinental  Ballistic  Missile  (ICBM)  booster  launched  the  first  Gemini  mission 
from  Cape  Canaveral  [4] .  This  program  led  to  the  first  American  space  walk;  and  in 
December  of  1965,  the  revolutionary  rendezvous  and  docking  concept  was  demon¬ 
strated  by  the  joint  rendezvous  mission  of  the  Gemini  VI  and  Gemini  VII  spacecrafts 
[5,  4],  Prior  to  the  conclusion  of  the  Gemini  program  in  1966,  work  had  already  be¬ 
gun  on  the  Apollo  program.  This  manned  program  used  the  rendezvous  techniques 
learned  during  the  Gemini  program  as  the  foundation  for  the  mission  to  the  Moon. 
In  the  Apollo  Moon  mission,  after  leaving  the  Lunar  surface,  the  Lunar  Excursion 
Module  (LEM)  was  required  to  rendezvous  and  dock  with  the  Command  Service 
Module  (CSM)  using  radar  during  long-range  separations  and  crew  observations  for 
the  final  docking  stages  [2],  Zimpfer  et  al.  [2]  explains  further  that  the  rendezvous 
calculations  were  performed  on  the  ground  while  the  on-board  system  was  able  to 
automatically  control  thruster  firings.  Again,  the  terminal  phase  of  capture  between 
the  LEM  and  CEM  was  manually  controlled  by  the  on-board  crew.  The  Apollo 
program  demonstrated  the  fundamentals  of  guidance,  navigation,  and  rendezvous  & 
capture,  with  minimal  computer  power  that  still  serves  as  the  foundation  for  today’s 
concepts. 
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1.1.2  Space  Shuttle  Program.  Perhaps  one  of  the  most  well  known  space 
vehicles  today,  the  Space  Shuttle  has  repeatedly  demonstrated  the  ability  to  ren¬ 
dezvous  in  LEO  and  perform  what  could  be  considered  the  start  of  proximity  opera¬ 
tions  since  the  1980’s.  Major  rendezvous  and  proximity  operations  accomplished  by 
the  Space  Shuttle  Program  include  the  Discovery  STS-51I  mission  in  1985,  where  the 
Shuttle  crew  captured,  repaired,  and  redeployed  the  SYNCOM  IV-4  (LEASAT-4) 
communications  satellite  [6].  In  years  following,  the  Orbiter  has  docked  with  the 
Soviet-made  Mir  space  station  and  the  International  Space  Station  (ISS).  In  gen¬ 
eral,  the  rendezvous  maneuver  sequence  is  planned  by  ground  operators,  and  with 
crew  command,  executed  by  the  on-board  Guidance  Navigation  &  Control  (GN&C) 
system  [2],  The  GN&C  system  performs  rendezvous  functions,  relative  navigation, 
targeting,  and  control;  yet,  the  Orbiter  crew  manually  performs  the  final  docking 
maneuvers  using  visual  aids  and  data  from  a  LIDAR  Trajectory  Control  Sensor  [2]. 
Including  the  Space  Shuttle  program  and  prior  programs,  the  level  of  computer  and 
sensor  technology  involved  in  actually  performing  the  rendezvous  and  proximity  op¬ 
erations  has  increased  dramatically.  Yet  the  final  stage,  the  action  defining  a  space 
proximity  mission,  still  remains  in  the  control  of  a  human  operator  performing  the 
final  action  using  today’s  available  computing  and  sensor  technology. 

1.1.3  Recent  Programs.  In  recent  years  there  have  been  several  programs 
designed  to  explore  autonomous  navigation  and  autonomous  proximity  operations 
and  planning.  Some  of  these  missions  include  the  Experimental  Satellite  System- 11 
(XSS-11)  of  the  USAF  and  the  Japanese  Hayabusa  mission  [2,  3]. 

The  USAF  Air  Force  Research  Laboratory  (AFRL)  and  Lockheed  Martin  Space 
Systems  Co.,  Waterton,  CO,  developed  the  XSS-11  spacecraft  to  demonstrate  au¬ 
tonomous  rendezvous,  autonomous  proximity  operations,  and  autonomous  mission 
planning.  Launched  by  a  Minotaur  I  booster  in  April  2005  from  Vandenberg  Air 
Force  Base,  the  spacecraft  successfully  performed  proximity  operations  with  its  ex¬ 
pended  upperstage  from  distances  between  500  m  and  1.5  km  [7].  The  spacecraft  has 
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performed  over  75  natural  motion  circumnavigations  of  the  Minotaur  upperstage  [7]. 
XSS-11  uses  advanced  on-board  planning  and  guidance  algorithms  including  an  event 
planner,  a  monitor  &  forward  thinking  resource  manager,  and  a  proximity  operations 
guidance  system  developed  by  Draper  Laboratory  [2],  The  100-kilogram  spacecraft 
is  equipped  with  an  active  LIDAR  sensor  and  a  passive  camera  used  for  relative  navi¬ 
gation  [2,  7].  Figure  1.1  shows  an  image  taken  by  the  30-kilogram  XSS-10  spacecraft 
(the  predicessor  of  the  XSS-11)  of  its  Boeing  Delta  II  upperstage  in  2003.  Figure  1.2 
is  from  an  image  taken  by  the  XSS-11  spacecraft  of  its  Minotaur  upperstage  from  a 
distance  of  500  meters  with  the  Earth  in  the  background  [7]. 


Figure  1.1  XSS-10  Image  of  Delta  II  Upperstage  (USAF  Photo) 


Figure  1.2  XSS-11  Image  of  Minotaur  Upperstage  (USAF  Photo) 
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In  May  9,  2003,  the  Hayabusa  (MUSES-C)  spacecraft  was  launched  with  a 
M-V  booster  from  the  Kagoshima  launch  site  by  the  Japan  Aerospace  Exploration 
Agency  (JAXA)  [8].  The  Hayabusa  spacecraft  was  developed  for  a  sample  return 
mission  to  the  asteroid  Itokawa  (25143),  named  after  the  late  Dr.  Hideo  Itokawa, 
considered  the  father  of  the  Japanese  space  program  [8].  The  mission  is  being  used 
to  demonstrate  technology  needed  for  the  spacecraft  to  autonomously  retrieve  sci¬ 
entific  samples  from  a  small  celestial  body  and  return  the  sample  back  to  Earth. 
JAXA  states  that  the  Hayabusa  mission  objectives  are  to  demonstrate  the  following 
four  key  technologies:  interplanetary  travel  using  ion  engines,  autonomous  naviga¬ 
tion  and  guidance  using  optical  measurements,  sample  collection  under  micro  gravity 
from  Itokawa ,  and  direct  re-entry  to  Earth  from  an  interplanetary  orbit  [3].  After 
successfully  navigating  itself  to  the  asteroid  without  human  guidance,  the  spacecraft 
began  its  proximity  operations  and  had  autonomously  observed  the  asteroid  by  hov¬ 
ering  7  to  20  kilometers  above  its  surface  [8].  Using  a  Telescope  Wide-View  Camera, 
LIDAR,  and  a  Near-Infrared  Spectrometer,  Hayabusa  has  mapped  the  surface  of 
the  asteroid  and  its  features  so  that  its  geometry,  mass  properties,  and  composition 
can  be  determined  [3,  8].  Figure  1.3  shows  an  image  of  the  asteroid  taken  by  the 
Hayabusa  spacecraft  camera  from  a  distance  of  approximately  8  kilometers  [3] . 


Figure  1.3  Asteroid  Itokawa  (courtesy  of  JAXA) 
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1.2  Applications 

In  a  paper  written  by  D.J.  Scheeres  [1],  he  explains  the  difficulty  of  a  spacecraft 
or  surface  vehicle  operating  in  the  vicinity  of  a  small  solar  system  body.  Mainly,  the 
difficulty  is  due  to  a  lack  of  knowledge  and  uncertainty  of  the  physical  characteris¬ 
tics  of  the  Target  body  and  its  possible  chaotic  dynamical  environment  [1].  Scheeres 
writes,  “To  successfully  carry  out  close  proximity  operations  about  these  bodies 
requires  an  understanding  of  the  orbital  dynamics  close  to  them,  a  knowledge  of 
the  physical  properties  of  the  body  and  the  spacecraft,  and  an  appropriate  level 
of  technological  sensing  and  control  capability  onboard  the  spacecraft  [1].”  Having 
the  ability  to  derive  and  estimate  physical  and  dynamical  properties  of  an  RSO 
based  on  post-rendezvous  observations  has  obvious  civilian  and  military  applica¬ 
tions.  In  all  cases,  the  pre-programmed  model  of  the  operating  environment  formed 
by  Earth-based  sensor  data  (including  the  dynamics  of  the  Target)  can  be  updated 
by  the  Observer  with  a  more  accurate  data  set  in  real-time;  thereby,  reducing  the 
uncertainty  of  the  operating  environment.  As  with  the  Hayabusa  mission,  a  civilian 
application  requiring  the  need  to  understand  the  physical  properties  and  dynam¬ 
ics  of  an  asteroid  or  small  moon  is  a  sample  return  mission.  Closer  to  Earth,  the 
inspection  of  space-based  assets,  such  as  the  ISS,  Hubble  telescope,  etc.,  could  be 
made  allowing  for  repair  or  servicing.  Possible  non-aggressive  military  applications, 
for  a  maneuverable  and  autonomous  spacecraft  capable  of  collecting  observational 
data  and  making  physical  and  dynamical  estimates  of  a  Target  range  from  dam¬ 
age  inspection  of  larger  US  satellites,  on-orbit  repair,  re-fuelling,  and  possibly  RSO 
identification  and  characterization. 

1.3  Key  Terms 

This  section  will  define  and  elaborate  on  key  terms  and  concepts  referred  to 
in  this  research.  A  general  definition  for  proximity  operations  and  dynamical  & 
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physical  characterization  is  given  below  as  well  as  clarification  on  nomenclature  used 
in  this  thesis. 

1.3.1  Space  Proximity  Operations.  In  general,  spacecraft  proximity  op¬ 
erations  is  defined  in  this  research  as  one  or  more  spacecraft  bodies  conducting 
an  activity  in  the  close  vicinity  of  a  RSO.  The  solution  to  the  relative  motion  of 
the  spacecraft  and  RSO  operating  in  close  proximity  to  each  other  is  given  by  the 
Clohessy-Wiltshire  equations  [9].  These  relative  motion  equations  along  with  their 
assumptions  are  explained  in  Chapter  2.  Although  no  specific  and  agreeable  defini¬ 
tion  was  found  in  the  literature  quantifying  the  terms  proximity  or  vicinity ,  dynam¬ 
ically,  Equation  2.5  must  be  satisfied  for  the  Clohessy-Wiltshire  equations  to  hold 
true.  It  would  be  acceptable  to  say  that  relative  distances  measured  in  the  order  of 
~  103  meters  or  less  could  be  classified  as  proximity  operations.  Traditionally,  the 
passive  body  of  interest,  in  which  the  observations  or  operations  are  being  conducted 
on,  is  termed  the  Target.  The  body  (i.e.  spacecraft)  that  is  actively  performing  the 
proximity  operations  (e.g.  observations)  is  termed  the  Chaser.  Other  terms  typi¬ 
cally  found  in  literature  refer  to  these  bodies  as  Chief  &  Deputy ,  Master  &  Slave , 
and  Target  &  Observer,  respectively.  This  paper  will  use  the  terms  of  Target  &  Ob¬ 
server;  where  the  Observer  is  a  spacecraft  performing  the  proximity  operation  and 
the  Target  is  the  Resident  Space  Object  of  interest.  Post-rendezvous,  the  Observer 
spacecraft  will  position  itself  in  the  vicinity  of  the  Target  in  preparation  to  carry 
out  its  proximity  operation  mission.  In  general,  proximity  operations  may  include 
specific  activities  such  as  taking  measurement  observations  of  the  Target,  docking 
with  the  Target  to  make  repairs,  or  performing  a  scientific  return  sample  mission. 
In  this  thesis,  the  proximity  activity  is  taking  observable  range  measurements  of  the 
Target  as  explained  in  later  sections. 

1.3.2  Spacecraft  Proximity  Sensors.  A  spacecraft  performing  proximity 
operations  would  be  equipped  with  a  sensor  suite  used  to  perform  standard  attitude 
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control  and  relative  navigation  which  would  be  needed  to  travel  to  the  point  of  ren¬ 
dezvous  where  the  proximity  mission  would  take  place.  This  sensor  suite  commonly 
includes  some  of  the  following:  a  sun  sensor,  earth  sensor,  star  sensor,  ring  laser 
gyro,  magnetometer,  and  Global  Positioning  System  (GPS)  unit  [10].  The  use  of 
a  star  sensor  (i.e.  star  camera)  is  a  method  of  far-range  photogrammetry  where 
the  distance  setting  of  the  source  is  essentially  infinity  and  object  measurements 
are  made  from  photographic  images  or  line-of-sight  measurements.  The  use  of  close- 
range  photogrammetry,  where  the  distance  of  the  source  is  finite,  can  be  employed  to 
determine  position  and  attitude  of  an  object;  where  as,  far-range  photogrammetry 
can  only  be  used  to  estimate  attitude  [11].  Close-range  photogrammetry  was  used 
for  the  the  XSS-10,  XSS-11,  and  Hayabusa  missions.  Using  photogrammetry  and 
pattern-recognition  methods,  position,  attitude,  and  physical  feature  information 
can  be  derived  from  the  imaged  target,  passively. 

Another  essential  sensor  for  space  rendezvous  and  proximity  activities  is  the 
scanning  Light  Detection  and  Ranging  (LIDAR)  instrument.  Rendezvous  Laser  Vi¬ 
sion  (RELAVIS)  and  the  Rendezvous  Laser  System  (RLS)  scanning  LIDAR  tech¬ 
nology,  built  by  Optech  Inc.,  Ontario,  Canada,  was  used  on  the  XSS-11  mission 
as  the  active  rendezvous  instrument  and  similar  LIDAR  technology  is  planned  for 
use  on  the  Mars  Sample  Return  mission  (MSR)  mission  as  the  primary  instrument 
for  autonomous  rendezvous  [12],  An  example  of  LIDAR  specifications  are  shown  in 
Table  1.1  [12], 

Table  1.1  LIDAR  RELAVIS  Instrument  Specifications  (Optech  Inc.) 


Min  Range 

0.5  m 

Max  Range 

5000  m 

Mass 

15  kg 

FOV 

10  deg 

Sample  Rate 

5  Hz 

Range  Accuracy 

1  cm 

Bearing  Accuracy 

0.35  mrad 

Attitude  Accuracy 

1  deg 
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A  scanning  time-of-flight  LIDAR  sensor  uses  a  pivoting  mirror  to  steer  laser 
pulses  sent  and  then  measures  the  optical  time-of-flight  required  for  the  pulses  to  re¬ 
turn  and  illuminate  the  focal  plane  array  along  with  measuring  the  angle  of  reflection. 
This  type  of  LIDAR  sensor  provides  range  and  bearing  of  the  target  and  can  also 
be  used  to  render  a  three-dimensional  point-cloud  image  of  the  RSO.  The  Hayabusa 
spacecraft  combined  LIDAR  range  data  with  camera  imagery  to  accurately  estimate 
the  spacecraft’s  position  relative  to  the  asteroid  [3].  Figure  1.4  shows  a  data  cloud 
image  of  the  asteroid  Itokawa  created  from  the  LIDAR  range  measurements  taken 
by  the  spacecraft  and  plotted  against  a  model  developed  by  S.  Ostro  et  al.  using 
ground-based  radar  [3]. 


Z  (km) 


0.3 


02 


Figure  1.4  LIDAR  Data  Cloud  of  Asteroid  Itokawa  (courtesy  of  JAXA) 

This  3-D  shape  model  was  developed  while  the  spacecraft  made  its  early  ob¬ 
servations  on  the  asteroid’s  equatorial  plane.  As  a  result,  insufficient  data  models 
the  characteristic  neck  of  the  asteroid  in  the  polar  regions  as  seen  by  the  optical 
navigation  camera  image  shown  earlier  in  Figure  1.3  [3].  Additional  data  fusion  will 
eventually  resolve  this. 
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To  explore  the  meaningfulness  and  possible  limitations  of  minimal  data  types, 
this  research  will  assume  that  the  only  available  data  from  the  sensors  of  the  Observer 
spacecraft  are  range  measurements  of  the  RSO  during  a  finite  period  of  time  when 
the  proximity  operation  occurs. 

1.3.3  Physical  &  Dynamical  Properties.  For  this  research,  the  physical 
properties  of  the  target  body  are  specifically  limited  to  the  moment-of-inertia  com¬ 
ponents  of  the  RSO.  The  dynamical  properties  are  characterized  by  the  relative 
position  and  velocity  states  as  well  as  the  angular  velocity  and  Euler  orientation 
angle  states.  Scheeres  et  al.  [13]  explains  that  using  ground-based  radar  observa¬ 
tions,  Ostro  et  al.  were  able  to  create  a  model  of  the  asteroid  Itokawa  and  estimate 
its  rotational  states  from  Earth.  Based  on  the  polyhedral  shape  model,  Scheeres  et 
al.  were  able  to  determine  estimates  of  the  asteroid’s  physical  properties  such  as  its 
dimensions,  volume,  center-of-mass  location,  and  moments  of  inertia  as  a  function 
of  the  unknown  mass  [13].  By  using  these  estimates  as  a  priori  knowledge,  a  space¬ 
craft  could  perform  observations  post-rendezvous,  re-estimate  the  dynamical  states 
of  the  RSO  at  a  significantly  closer  range  (improving  observational  accuracy),  and 
update  the  pre-programmed  earth-based  data  set  resulting  in  a  high-fidelity  model 
and  accurate  state  estimates.  Physical  properties  such  as  mass,  moments  of  inertia, 
geometry  could  be  validated  and  well  characterized  to  relatively  high  accuracy.  With 
an  updated  high-fidelity  model  of  the  RSO,  and  with  new  rotational  state  informa¬ 
tion,  understanding  the  dynamical  properties  of  the  RSO  (e.g.  axis  of  rotation  and 
spin  stability)  is  greatly  enhanced. 

1-4  Literature  Review 

This  thesis  primarily  involves  spacecraft  dynamics  and  estimation  theory.  Lit¬ 
erature  is  abundant  for  both  relative  motion  dynamics  and  rotational  kinematics  for 
a  rigid-body.  Texts  by  W.  Wiesel  [9]  and  B.  Wie  [14]  cover  both  relative  motion  and 

rigid-body  dynamics  quite  well.  For  this  research,  those  references  are  the  primary 
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source  for  the  dynamics  theory  developed  in  Chapter  2.  Another  text  written  by  W. 
Wiesel  [15]  provides  a  very  complete  foundation  and  guide  for  the  implementation  of 
least  squares  estimation  theory  for  orbit  determination,  which  is  heavily  used  in  this 
research.  In  addition,  numerous  journals  authored  by  D.  Scheeres  [13,  1]  provides  in¬ 
sight  into  applications  of  estimation  techniques  using  LIDAR  and  photogrammetry, 
proximity  missions,  and  methods  in  characterizing  physical  and  dynamical  proper¬ 
ties  of  an  asteroid.  Finally,  theses  authored  by  T.  Storch,  J.  Heslin,  and  V.  Chioma 
utilized  least  squares  estimation  methods  for  various  orbit  estimation  applications. 
However,  methods  in  directly  estimating  the  moments  of  inertia  of  a  rigid  body  based 
solely  on  range  observations  of  the  rotating  body  was  not  found  in  this  literature 
review. 

1.5  Research  Objective 

This  thesis  will  explore  the  ability  of  a  nonlinear  least  squares  algorithm  to 
estimate  the  relative  motion,  rotational  states,  and  moments  of  inertia  of  a  RSO  by 
processing  range  data  collected  from  tracked  observational  points  on  the  RSO  body 
from  a  simulated  proximity  mission. 
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II.  Theoretical  Basis 


2. 1  Background 

This  chapter  will  discuss  the  theory  that  is  the  foundation  of  the  methodology 
used  in  this  thesis.  The  theoretical  basis  of  this  research  involves  relative  motion  and 
rotational  dynamics,  linearization  of  nonlinear  dynamical  systems,  and  estimation 
theory;  specifically,  the  nonlinear  least  squares  method. 

2.2  Dynamics 

This  section  defines  the  reference  frames  used  and  develops  the  equations  of 
motion  necessary  to  describe  relative  motion  and  rotational  kinematics  following 
derivations  provided  by  Wiesel  [9],  Wie  [14],  Meirovitch  [16],  and  Curtis  [17].  This 
section  results  in  the  formation  of  the  dynamics  state  matrix  for  the  Target  body  as 
seen  by  the  Observer  spacecraft. 

2.2.1  Reference  Frames.  The  coordinate  frames  used  in  this  research  are 
all  Cartesian  dextral  frames.  Along  with  each  frame  is  a  set  of  basis  vectors  and 
their  associated  scaler  coordinates. 

2. 2. 1.1  Earth-Centered  Inertial.  The  Earth-Centered  Inertial  (ECI) 
reference  frame  is  denoted  as  31:  with  the  origin  at  the  center  of  mass  (CM)  of  the 
Earth.  The  31  frame  has  basis  vectors  of  {I,  J,  K}  where,  the  I  unit  vector  points  in 
the  direction  of  the  vernal  equinox,  traditionally  denoted  by  7.  The  K-axis  is  defined 
as  the  axis  aligned  with  the  spin-axis  of  the  Earth  and  points  out  of  the  geographic 
north  pole.  The  J-axis  completes  this  orthogonal  frame  shown  below  in  Figure  2.1. 

The  scalar  coordinates  associated  with  the  basis  vectors  of  this  frame  are 
( X,Y,Z ).  The  3 i  frame  is  considered  inertial  and  is  therefore,  fixed  in  space  and 
does  not  rotate  with  the  Earth. 
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2.2. 1.2  Orbital  Frame.  The  orbital  reference  frame,  denoted  as 
is  also  referred  to  as  the  Euler-Hill  frame  in  most  relative  motion  literature.  The 
rotating  orbital  frame  has  basis  vectors  of  {§r,  eg,  ey}  where  the  origin  of  lies 
on  a  circular  orbit  defined  by  the  vector,  R tgt,  measured  from  the  &  origin  at  the 
center-of-mass  of  the  Earth.  The  orbital  frame  rotates  about  the  ejy-axis  with  the 
mean  orbital  motion,  n,  defined  later  in  this  chapter.  The  e#-axis,  referred  to  as  the 
in-track  direction,  points  along  the  velocity  vector  of  the  motion  and  is  tangent  to 
the  circular  orbit  during  the  entire  orbital  period.  The  ey-axis  is  perpendicular  with 
the  orbital  plane  and  is  termed  the  cross-track  direction.  The  e«-axis  points  radially 
outward,  in  the  radial-track  direction,  the  same  direction  as  and  completes  the 
dextral  reference  frame  shown  in  Figure  2.2. 

The  center  of  mass  of  the  Target  body  is  centered  with  the  origin  of  There¬ 
fore,  the  vector  Rtgt  points  from  the  origin  of  5,:  to  the  origin  of  5,  -where  the  CM  of 
the  Target  is  located.  It  is  important  to  note  that  the  Target  body  has  freedom  of 
rotation,  with  the  constraint,  that  the  CM  of  the  Target  remains  fixed  on  the  $e  ori¬ 
gin.  The  scalar  coordinates  associated  with  the  basis  vectors  of  $e  are  (8x,8y,8z). 
Note  that  Sx  and  Sy  do  not  necessarily  have  to  represent  rectilinear  coordinates. 
They  can  be  interpreted  as  cylindrical  curvilinear  coordinates  where  Sx  represents 
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Figure  2.2  Orbital  Frame  of  Target,  $e 

a  radial  displacement,  Sr,  and  8y  represents  the  curved  flight-path  displacement, 
r0bs59 ,  without  change  to  the  mathematical  expressions  [18,  9].  Using  a  curvilin¬ 
ear  coordinate  system  does  not  restrict  the  magnitude  of  the  e#  displacement;  and 
therefore,  extends  the  validity  of  the  relative  motion  solution  [9,  14].  However,  in 
this  research  the  Target-Observer  distances  are  so  small  that  differences  between  the 
two  frames  are  considered  negligible;  therefore,  no  distinction  will  be  made  between 
the  linear  and  curvilinear  frames. 

2. 2. 1.3  Body  Frame.  The  body  reference  frame  of  the  Target,  de¬ 
noted  as  is  non- inertial  and  fixed  to  the  body  of  the  Target.  The  origin  of  $7, 
is  placed  at  the  center-of-mass  of  the  Target.  The  frame  has  basis  vectors  of 
{bl5b2,b3}  which  are  aligned  with  the  principal  axes  of  inertia  of  the  Target  body. 
Figure  2.3  illustrates  a  generic  RSO  and  its  body  frame.  The  associated  scalar 
coordinates  are  (. bi,b2,b3 ). 

2.2.2  Relative  Motion  Dynamics.  Originally  introduced  by  Hill  in  1878, 
W.  H.  Clohessy  and  R.  S.  Wiltshire  developed  a  linearized  set  of  equations  of  relative 
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Figure  2.3  Body  Frame  of  Target, 

motion  in  1960  [19].  Termed  the  Clohessy-Wiltshire  (CW)  equations,  also  known  as 
Hill’s  equations,  they  describe  the  relative  motion  of  two  satellites  in  close  proximity 
while  one  satellite  is  in  a  circular  orbit  around  a  central  body.  The  classic  CW 
equations  that  are  derived  in  this  section  assume  the  Target  body  follows  a  circular 
reference  orbit  around  the  Earth  and  neglects  any  disturbance  forces.  It  should  be 
mentioned,  that  the  CW  equations  have  been  modified  by  others  for  application  in 
eccentric  reference  orbits  and  to  also  include  perturbations. 


2.2.2. 1  Hill- Clohessy-Wiltshire  Equations.  The  governing  equation 
of  motion  for  the  classic  two-body  problem  of  a  satellite  and  central  body  is 


r  = 


(2.1) 


where  r  is  the  acceleration  of  the  body  and  r  is  the  position  vector  of  the  body. 
Applying  two-body  motion  to  the  Target  in  circular  orbit  around  the  Earth  results 
in  Equation  2.1  re-written  as 


R 


■tgt 


/^©R  tgt 
—R 3 

1Xtgt 


(2.2) 


where  /%  is  the  gravitational  parameter  of  the  Earth. 
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The  position  vector  of  the  Target  Rtst  from  the  origin  of  5,  is  dehned  as 


R-tgt  —  Rtgt&R  (2-3) 

and  in  conducting  a  space  proximity  operation,  consider  the  Observer  spacecraft 
operating  in  close  vicinity  to  the  Target.  Keep  in  mind  that  the  center  of  mass 
of  the  Target  is  fixed  to  the  origin  of  the  rotating  orbital  frame,  $e.  The  relative 
position  vector  of  the  Observer  from  the  Target,  8p,  is  defined  by  the  equation  below 

Sp  =  8xeR  +  8yeg  +  8zeN  (2.4) 


A  major  assumption  is  that  the  magnitude  of  dp  is  relatively  small  compared  to  the 
magnitude  of  Rf(/t .  That  is 


M 

R-tfjt  | 


<C  1 


(2.5) 


With  that  assumption  in  mind,  the  Observer  follows  two-body  motion  as  well 


Yobs 


r 


3 

obs 


(2.6) 


where  the  position  vector  from  the  origin  of  to  the  Observer  spacecraft  is  r0i)S . 
The  vector  relation  between  the  Observer  and  Target  is  then  given  by  the  relation 
below  and  shown  by  the  exaggerated  illustration  in  Figure  2.4. 


obs  R  tgt  +  8  p 


(2.7) 


Substituting  Equations  2.3  and  2.4  in  Equation  2.7  results  in 

^ obs  ( Rtgt  +  8x)eR  8yee  8z&]y  (2.8) 
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X  / 


Figure  2.4  Relative  Motion  System 

The  following  acceleration  formula  [9,  17]  along  with  the  circular  orbit  assumption 
r*  =  u)ei  x  re  +  ujet  x  (u)ei  x  re)  +  2 'wei  x  re  +  re  (2.9) 


is  used  where  the  superscripts  i  and  e  indicate  the  reference  frames  and  $e, 
respectively.  The  inertial  acceleration  of  the  Observer  spacecraft  is  found  by  taking 
the  first  and  second  derivatives  of  Equation  2.8  with  respect  to  time  and  making  the 
proper  substitutions  into  Equation  2.9  resulting  in 

robs  =  [5x  -  2t o5y  -  (jJ2(Rtgt  +  +  ( Sy  +  2u5x  -  u25y)ee  +  SzeN  (2.10) 

where  u  is  the  angular  velocity  between  the  frames  $e  and  With  the  assumption 
that  the  Target  body  follows  a  circular  orbit,  the  orbital  frame  rotation  is  equal  to 
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the  mean  orbital  motion  of  the  Target 


ojei  =  Lue^  =  ne^  or  simply  u  =  n 


(2.11) 


where  the  mean  orbital  motion  of  the  Target,  n,  is  defined  as 


n  = 


Me 

i?3 

tgt 


(2.12) 


Substituting  Equation  2.8  into  Equation  2.6,  now  results  in 


f* obs 


y^Rtgt  +  Sx)eR  +  Syee  +  SzeN\ 
[(Rtgt  +  Sx)2  +  5y2  +  5z2]  § 


(2.13) 


Using  the  binomial  theorem  in  the  form  of 


(1  +  b)n  =  1  +  nb  + 


n(n-l)  2 

2! 


(2.14) 


the  denominator  of  Equation  2.13  is  expanded  as 


r“3  = 

'  obs 


[(^? tgt  2Rtgtfix  +  fix  +  Sy  +  Sz  ]  2 

Sx 2  5y2  hz2 


Rtgt  (  1 


-5a; 


_  p— 3 

— 


Rtgt 

3/2  <hr 

n 

tgt 


R2 

11 tgt 


H - - - b 

^  /?2  ^ 
tgt 


Rtgt  /  J 


2  V#; 


5a;2  5y2  hz2 

+  — 6  — 6  — I-  ••• 


Kt 


K. 


(2.15) 

(2.16) 

(2.17) 


Replace  the  denominator  of  Equation  2.13  with  the  binomial  expansion  and  neglect 
terms  higher  than  first-order.  This  results  in 


obs  ~  733  \.(Rtgt  +  SyGg  ~\~  SzQjy] 

Ktgt 


(2.18) 
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Equate  Equation  2.10  with  Equation  2.18,  substitute  Equations  2.11  and  2.12,  re¬ 
sulting  in  the  system  of  equations 


Sx  —  2  ujSy  —  u2(Rtgt  +  Sx) 

Rfgt  2$X 

0 

Sy  +  2u>Sx  —  u)2Sy 

+  n2 

Sy 

= 

0 

Sz 

Sz 

0 

Solve  each  equation  (i.e.  row)  of  Equation  2.19  for  Sx,  Sy,  and  Sz,  respectively. 
Finally,  the  system  of  CW  equations  [9]  describing  the  relative  motion  of  the  Observer 
spacecraft  with  respect  to  the  Target  is  given  by 

Sx  =  3  n2Sx  +  2  nSy  (2.20) 

Sy  =  —2  nScc  (2-21) 

Sz  =  —n2Sz  (2.22) 


to  the  CW 
Sx(t ) 

Sz(t) 


2. 2. 2. 2  Hill-  Clohessy-  Wiltshire  Solution.  The  closed-form  solution 
Equations  [9,  17]  is  given  below 


=  (4  —  3  cos  nt)Sx  o  + 


sin  nt 


n 


=  6(sin  nt  —  nt)5x  0  +  Sy0  — 

=  fc0coSnt  +  h£smni 
n 


Sx o  -I — (1  —  cos  nt)  Syo 
n 

2  4sinnt  —  3nt 

—  (1  —  cos  nt)dx0  H - 

n  n 


(2.23) 

^o(2.24) 

(2.25) 


Recall  that  the  relative  position  vector  of  the  Observer  spaceraft  with  respect  to  the 
Target  in  the  $e  frame  is  Sp  (Equ.  2.4).  Therefore,  the  relative  velocity  vector,  Sp, 
is  computed  by  taking  the  first  time-derivative  of  the  position  vector;  defined  as 


Sp  = 


Sx 

/  \ 

Sx 

Sy 

>  and  Sp  =  < 

Sy  > 

Sz 

) 

Sz 

\  / 

(2.26) 
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Therefore,  the  analytical  solution  for  the  relative  velocity  components  are  found  by 
taking  the  time-derivative  of  Equations  2.23,  2.24,  and  2.25,  resulting  in  the  following 
three  expressions: 


5x(t)  = 

Sx o  cos  nt  +  (3n6xo  +  2<b/o)  sin  nt 

(2.27) 

II 

(6nSx0  +  4:Sy0)  cos  nt  —  2 Sx0  sin  nt  —  6 n5x0  —  38y0 

(2.28) 

II 

—5zon  sin  nt  +  Szq  cos  nt 

(2.29) 

The  CW  state  vector  explicitly  written  as  a  function  of  time  is  then  defined  as 


x»(t) 


Sx(t ) 
8y(t) 
5z(t ) 
5x(t ) 
8y(t) 

5z(t ) 


(2.30) 


Let  Xcu,(to)  represent  the  initial  CW  state  at  an  initial  epoch  time,  to.  Then  let 
Xcw(t)  represent  the  propagated  CW  states  at  some  future  time,  t.  Using  both  the 
closed-form  CW  solution  for  the  relative  position  and  velocity,  define  the  CW  state 
transition  matrix,  as  the  matrix  to  propagate  from  XCU)(t0)  to  Xcw(t),  as 


^ cw (U  f'o) 


4  —  3  cos  nt 
6  (sin  nt  —  nt) 

0 

3 n  sin  nt 
— 6n(l  —  cos  nt) 
0 


0  0 

1  0 

0  cos  nt 

0  0 

0  0 

0  —  n  sin  nt 


-sin  nt  -(1  —  cosnt)0 

— 1(1  —  cos  nt)  ^(4sinnt  —  3nt)  0 

0  0  -sin  nt 

n 

cos  nt  2  sin  nt  0 

—2  sin  nt  4  cos  nt  —  3  0 

0  0  cos  nt 


(2.31) 
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Finally,  the  complete  solution  of  the  relative  motion  system  in  the  frame  is  com¬ 
pactly  defined  as 

Xcu,(t)  =  ^cw(t,t0)Xcw(t0)  (2.32) 

2.2.3  Rotational  Dynamics.  The  rotational  equations  of  motion  of  a  rigid 
body  are  given  by  two  systems  of  equations.  The  first  system  is  a  set  of  three 
differential  equations  known  as  Euler’s  Equations.  The  second  set  consists  of  three 
differential  equations  describing  the  orientation  of  the  body  in  terms  of  the  classic 
Euler  angles.  It  should  be  mentioned  that  a  singularity  does  exist  when  using  Euler 
angles  (which  will  be  discussed  later  in  this  chapter).  One  could  also  obtain  the 
second  set  of  equations  of  motion  using  quaternions  to  describe  the  rigid  body’s 
orientation  (without  singularities).  However,  for  the  scope  of  this  research  it  is 
assumed  that  no  operation  is  taking  place  at  the  singularity.  It  is  also  assumed  that 
the  Target  body  is  a  rigid  body. 

2.2.3. 1  Euler’s  Rotational  Equations  of  Motion.  Euler’s  Equations, 
as  derived  in  [9]  and  [14]  begins  with  the  relation  that  the  applied  torque  must  equal 
the  inertial  time-rate-of-change  of  the  angular  momentum  of  a  rigid  body  (i.e  the 
Target  body) 

M  =  H  (2.33) 

where  M  is  the  external  moment  acting  about  the  CM  of  the  Target  body  and  H  is 
the  angular  momentum  also  about  the  CM  of  the  Target. 

The  angular  momentum  of  the  Target,  in  the  body  frame,  is  expressed  as 

H  =  Iojbi  (2.34) 

Since  this  is  a  rigid  body,  the  moment  of  inertia  matrix,  I,  is  constant  in  Then 
using  what  is  sometimes  referred  to  as  the  transport  theorem,  the  time  derivative  of 
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the  angular  momentum  in  the  inertial  frame  is  given  by 


H  = 


+  u>bi  x  H 


(2.35) 


where  the  b  and  i  superscripts  refer  to  the  and  &,  respectively,  where  the  time- 
derivative  or  quantity  is  in  reference  to.  The  angular  velocity,  u>b\  is  simply  the 
angular  velocity  of  with  respect  to  the  and  it  is  also  the  same  quantity  used 
in  Equation  2.34.  Substituting  Equation  2.34  into  Equation  2.35  and  recalling  that 
the  moment-of-inertia  matrix  of  the  body  is  constant  in  57  yields  Euler’s  rotational 
equation  of  motion  in  vector  form 


M  =  Iujbi  +  ojbi  x  Iojbi 


(2.36) 


The  term  Cjbl  is  time-derivative  of  the  angular  velocity  in  the  body  frame. 

Since  the  57  frame  is  defined  as  the  principle- axis  reference  frame  of  the  Target,  the 
moment-of-inertia  matrix  of  the  Target  body  is  given  by 


I  = 


A  0  0 
0  5  0 
0  0  C 


(2.37) 


in  the  5&  frame.  The  angular  velocity  vector,  u>bl  is  defined  by 


UJbi  —  uqbi  +  (^2^2  +  UTjbs 


(2.38) 


In  this  research,  no  external  torques  are  applied  to  the  Target.  Therefore, 

M  =  Obx  +  0b2  +  0b3  (2.39) 
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Substituting  Equations  2.37  through  2.39  into  Equation  2.36  in  component  form 
yields  Euler’s  Equations,  as  derived  in  [9],  [14],  and  [16],  shown  below 


0  =  A(jJi  ( C  — 

(2.40) 

0  =  Bd)2  ~\~  {A  — 

(2.41) 

0  =  Cd)  3  +  (yB  —  jA^(jJ\U02 

(2.42) 

Solving  for  the  angular  accelerations  results  in  the  following  coupled,  nonlinear,  first 

order  differential  equations  of  motion 

*1  =  (B~C)LU2La3 

(2.43) 

( C-A ) 

(2.44) 

(- A-B ) 

w, 3  —  i0\i02 

(2.45) 

2. 2. 3. 2  Rotational  Kinematics.  Completing  the  set  of  equations  of 
motion  for  a  rigid  body  are  the  kinematic  differential  equations  describing  the  orien¬ 
tation  of  the  Target.  A  3-2-1  Euler  angle  sequence  (i.e.  roll-pitch-yaw)  is  used  and 
derived  following  Wiesel  [9],  Wie  [14],  and  Meirovitch  [16]. 

The  3-2-1  Euler  sequence  describes  the  rotation  sequence  of  a  body  frame 
initially  coincident  with  an  inertial  reference  frame  rotated  to  a  final  orientation 
described  by  the  angles,  ffi,  92,  63.  In  the  case  of  this  research,  the  rotating  frame 
is  considered  quasi-inertial  with  respect  to  the  $i,  frame  for  short  and  finite  periods 
of  time  during  the  proximity  operation  and  acts  effectively  as  a  reference  frame  for 
the  orientation  of  the  Target  body  frame. 

Beginning  with  the  standard  sequence  to  rotate  from  the  orbital  reference  frame 
to  the  body  frame,  the  matrix  [R]  is  a  rotation  matrix  of  sines  and  cosines.  The 
subscript  of  [R]  indicates  the  axis  in  which  a  rotation  of  the  specified  angle  will 
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occur.  The  basis  vectors  {£,?y,C}  and  {£/,t/>C/}  describe  the  intermediate  frames 
during  the  rotation  sequence. 


L  [R3(03)]  :  {eR,e0,eN}  ->  {£',  v'X'} 

2.  \R2(d2)]:{?,rf,C'}^{Z,ri,C} 

3.  [Ri^r)]  :{e,»7,C}-^{bi,b2,b3} 

Mathematically,  this  transition  between  frames  occurs  as  such 


j  v  — 


=  [R3(03)K 


=  F*(02)K 


=  [Ri(^i)]U 
C 


where  the  elementary  matrices  are  defined  as 


PM»s)] 


[R2(#2)] 


[Ri(»i)] 


cos  03  sin  03  0 

—  sin  03  cos  03  0 

0  0  1 

cos  6*2  0  —  sin  62 

0  1  0 

sin  62  0  cos  O2 

1  0  0 

0  cos  0i  sin  0i 

0  —  sin  0i  cos  0i 
2-13 


(2.46a) 


(2.46b) 


(2.46c) 


(2.47a) 


(2.47b) 


(2.47c) 


Equating  the  intermediate  transformation  equations,  Equations  2.46a-c,  yields  the 
relation 


/  \ 

bi 

/  \ 

&R 

b2 

>  =  [Ri(di)]  [R2(02)]  [R3(03)]  < 

ee  > 

b3 

V  > 

\  / 

(2.48) 


By  performing  matrix  multiplication  in  the  proper  sequence  for  each  elementary 
rotation  matrix  yields  a  single  rotation  matrix  denoted  as  [R6e] 


[Rr (0i)]  [R2(02)]  [R3(03)] 

c92c  6*3  c92s93  -s02 

s6is62c03  —  c$is#3  s6\s92s9^  +  c9\c9s  s0ic02 
c9is92c9s  +  s#isd3  c9is92s93  —  s$icd3  c9ic92 


(2.49) 

(2.50) 


where  the  abbreviations  s 9  =  sin($),  and  c 9  =  cos (0)  are  used.  The  superscript 
above  the  rotation  matrix  be  denotes  transformation  from  — >  $i,.  Therefore,  a 
vector  described  in  the  orbital  frame,  denoted  as  can  be  transformed  into  a 

vector  written  in  body  frame  components,  {^},  by  pre- multiplying  { $e }  with  [R6e] 
as  shown 

{&}  =  [R'l  {&}  (2.51) 

The  rotation  matrix,  [R6e] ,  is  an  orthonormal  matrix;  hence,  it  can  be  shown  that 
the  inverse  of  the  rotation  matrix  is  simply  its  transpose.  Therefore,  the  rotation 
matrix,  [Refe],  for  the  inverse  transformation,  ^6  — * ►  is  simply 

[Keb]  =  [R6e]T  (2.52) 


Since  the  focus  of  the  this  paper  is  based  on  the  relative  motion  of  the  Target  and 
the  Observer  with  respect  to  ,  it  will  be  necessary  to  describe  the  orientation 
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of  the  Target  body  with  respect  to  the  orbital  frame  as  well.  Therefore,  a  vector 
written  in  the  body  frame,  {St,},  is  transformed  into  a  vector  written  in  orbital  frame 
components,  {3>  },  using  the  rotation  matrix,  [Re6] ,  as  shown 

{&}  =  [Re6]  {&}  (2.53) 

where 

C$2  C$3  s9iS92c93  —  C$1  S#3  c9iS92c93  +  sdisd3 

[Re6]  =  c92s93  s9is92s93  +  c93c93  c9is92s93  —  s93c93  (2-54) 

—s92  s$ic#2  c$i cd2 

To  obtain  the  kinematic  equations  of  motion,  the  time-dependence  of  the  Euler  angles 
will  be  formed.  The  time  derivatives  of  the  3-2-1  Euler  angles,  termed  Euler  rates, 
are  denoted  as  9 92,  and  93. 

The  angular  velocity  solution  of  the  body  is  provided  by  numerical  integration  of 
Euler’s  Equations  2.43,  2.44,  and  2.45.  The  angular  velocity  vector  is  equal  to  the 
angular  velocity  of  the  body  frame  with  respect  to  the  inertial  frame  as  shown  earlier 
in  Equation  2.38. 

u>  =  ujbe  (2.55) 

Therefore,  with  the  rotational  sequence  in  mind,  one  can  write  the  equivalent  state¬ 
ment  with  mixed  basis  vectors 

u}  =  d]bi  +  92r\  +  93Q  (2.56) 

This  relation  describes  the  angular  rate  of  rotation  for  each  axis  during  the  coordinate 
frame  rotation  sequence.  Using  Equation  2.56  and  Equations  2.47a-c  results  in  the 
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angular  velocity  expressed  solely  in  the  body  frame 


/  \ 

CJi 

/  \ 

M 

/  \ 

0 

<  Cc>2 

>  =  < 

0 

>  +  [1^1  (^l)]  < 

e2  \  +  [Ri(0r)]  [r2(02)]  < 

0  > 

UJ3 

V  / 

0 

k  / 

1*>J 

03 

V  / 

(2*57) 


The  result  of  the  multiplication  in  matrix  form  is 


/  \ 

u  1 

10  —  sin  d2 

/  \ 

0i 

< 

>  = 

0  cos  61  sin  6\  cos  02 

< 

02  > 

LV3 

V  > 

0  —  sin  6\  cos  6*i  cos  02 

03 

V  / 

(2.58) 


Finally,  solving  for  the  Euler  rates  yields  the  remaining  system  of  kinematic  dif¬ 
ferential  equations  [14]  needed  for  describing  the  rotational  motion  of  the  Target 
body 

1 

COS  02 

and  expressed  in  three  scalar  equations  for  consistency 


COS  62 

sin  6*i  sin  d2 

cos  6\  sin  62 

CJi 

0 

cos  6*i  cos  02 

—  sin  61  cos  6*2 

< 

U2  * 

0 

sin  61 

cos  61 

U3 

(2.59) 


6\  =  ui  +  UO2  sin  6*i  tan  02  +  <^3  cos  6\  tan  02  (2.60) 

02  =  1-^2  cos  6\  —  sin  0\  (2-61) 

03  =  u2  sin  61  sec  02  +  uq  cos  61  sec  02  (2.62) 


The  solution  to  Euler’s  Equations  2.43,  2.44,  and  2.45,  provides  the  angular  velocity 
vector  as  a  function  of  time.  Using  that  vector  in  Equation  2.59,  and  numerically 
integrating,  will  provide  the  orientation  of  the  Target  as  a  function  of  time. 
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The  major  drawback  of  using  Euler  angles  to  describe  the  orientation  of  a  rigid  body 
is  the  existence  of  a  singularity.  In  this  rotational  scheme,  the  singularity  occurs 
when  62  =  | .  As  mentioned  in  the  beginning  of  this  section,  for  the  purpose  of  this 
research,  Euler  angles  are  sufficient  to  describe  the  orientation  of  the  Target.  It  will 
assumed  that  the  Target  body  does  not  operate  in  this  singularity  region  during  the 
relatively  short  duration  of  the  proximity  operation  studied. 


2.2.4  Dynamics  State  Vector.  The  Target  is  characterized  by  its  states, 
including  both  its  relative  motion  (i.e.  CW  position  and  velocity)  with  the  Observer 
as  well  as  its  rotational  body  dynamics  (i.e.  angular  rotation  and  Euler  orientation 
angles).  Consolidating  the  state  parameters  that  describe  the  dynamical  uniqueness 
of  the  Target  at  any  moment  in  time  during  the  proximity  operation  leads  to  the 
formation  of  the  complete  dynamics  state  matrix. 

Define  Xrot(f)  as  the  rotational  dynamics  portion  of  the  state  matrix 


Xrot(t)  — 


uq(f) 

w2(i) 

^3  (t) 

°l(t) 

e2(t) 


(2.63) 


Using  a  Runge-Kutta  numerical  integration  method,  the  numerical  solution  of  Equa¬ 
tions  2.43,  2.44,  and  2.45  results  in  uq,  cn2,  and  u>3  and  the  numerical  solution  of 
Equations  2.60,  2.61,  and  2.62  results  in  6\,  d2,  and  03.  These  parameters  of  angular 
velocity  and  orientation  angles  form  a  state  vector  of  the  rotational  states  of  the 
Target  as  a  function  of  time. 

Recall,  Equation  2.30  defines  the  states  for  the  dynamics  of  the  relative  motion 
between  the  Target  and  the  Observer.  Consolidating  the  CW  and  rotational  state 
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vectors  results  in  the  complete  state  vector,  X(t),  describing  the  dynamics  of  the 
Target  body  with  respect  to  time. 


X(t)  = 


XCU)(t) 

Xrot(t) 


(2.64) 


In  expanded  form,  the  transpose  of  the  complete  12-state  vector  of  the  Target’s 
dynamics  is  defined  as 


X  = 


Sx  Sy  Sz  Sx  Sy  Sz  uq  u>2  <^3  O2 


(2.65) 


Using  the  assumption  that  the  moments  of  inertia  of  the  Target  are  constant  over 
time.  The  time-derivative  of  the  Target  MOI  are  given  by  the  following  equations 


A  =  0 


B  =  0 


<7  =  0 


where  the  MOI  states  are  defined  as 


y 

moi 


(2.66) 

(2.67) 

(2.68) 


(2.69) 


By  removing  the  CW  states  XC)(!  from  Equation  2.64  and  including  X„t0( (t),  results 
in  a  separate  9-state  vector  version  defined  as 


X  = 


l 0\  ui2  ^3  9\  &2  O3  A  B  C 


(2.70) 
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2. 3  Dynamics  Linearization 

Apparently,  the  world  was  not  created  using  a  linear  dynamics  model,  which 
simply  means  most  of  the  dynamical  behavior  that  one  would  be  interested  in  is  sim¬ 
ply  nonlinear  and/or  unpredictable  in  nature.  However  by  limiting  our  observation 
time  and  invoking  small  state  changes,  one  can  essentially  linearize  the  nonlinear 
systems  in  regions  that  would  provide  approximate  results.  Using  the  derivations 
provided  in  two  texts  by  Wiesel  [15,  20],  methods  of  linearization  are  formulated 
below. 

2.3.1  Equations  of  Variation.  The  dynamical  state  of  the  Target-Observer 
system  is  described  by  the  relative  position,  relative  velocity,  angular  velocity  com¬ 
ponents,  and  Euler  orientation  angles  derived  in  the  last  section.  These  system  state 
variables  comprise  the  dynamics  state  vector,  X.  Each  system  state  variable  has  an 
associated  equation  of  motion.  Therefore,  the  change  of  the  state  vector  with  respect 
to  time  is  dictated  by  the  equations  of  motion 

X  =  f(X,t)  (2.71) 

To  create  a  trajectory,  initial  state  conditions,  X0(t0),  must  be  defined  at  some 
epoch  time,  t0.  Using  numerical  integration  methods,  the  state  vector  can  then 
be  propagated  to  some  future  time  resulting  in  a  trajectory  X0(£).  With  a  slight 
change  in  the  state  vector,  <5X,  one  would  expect  a  slight  change  in  the  trajectory. 
The  resulting  trajectory  is  then  said  to  be  close  to  the  original  trajectory.  This 
fundamental  expectation  can  be  written,  in  general,  as 

X(f)  =  Xo(t)  +  5X(t)  (2.72) 
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Substituting  Equation  2.72  into  Equation  2.71  yields 


X  =  X0  +  SX  =  f  (Xo  +  5X,  t ) 


(2.73) 


Applying  a  Taylor  series  expansion  to  the  equations  of  motion 

X  !—P~(x  ~  “)”  =  /<0|(a)  +  fm(a)(x  -  a)  +  , , ,  (2,74) 

71=0 

about  the  point  a  =  X0  +  AX  where,  AX=0,  and  neglecting  second-order  terms  and 
higher,  results  in 

X  «  f(X0  +  0)  +  f(1)(X0  +  0)(X0  +  SX  -  X0  -  0)  +  0(2)  (2.75) 


X0  +  SX  «  f  (X0)  + 


df 

dX 


5X 


Xo 


(2.76) 


After  cancelling  Equation  2.71  from  both  sides,  results  in  the  familiar  equations  of 
variation  used  for  the  linearization  of  the  dynamics 


5X  =  A(t)5X 

Xo 


(2.77) 


The  equations  of  variation  are  time-varying  linear  ordinary  differential  equations. 
Matrix  A  is  a  square  matrix  with  the  dimensions  of  the  number  of  states.  It  is 
formed  by  the  partial  derivatives  of  the  equations  of  motion  with  respect  to  the 
state  variables  evaluated  at  some  nominal  trajectory. 


2.3.2  State  Transition  Matrix.  The  solution  to  the  Equations  of  Variation, 
Equation  2.77,  is  provided  in  detail  by  Wiesel  [15].  They  are  found  by  summing  the 
individual  components  of  SX  at  an  initial  time  multiplied  by  each  solution  vector, 
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that  is 


N 

SX(t)  =  Y,SXi(t0)<f>i(t)  (2.78) 

1=1 

Arranging  the  vector  solutions,  </>, ,  into  a  square  matrix  format  results  in  the  forma¬ 
tion  of  the  State  Transition  Matrix ,  $ 

5X(t)  =  <f>(t,t0)5X(t0)  (2.79) 

The  matrix  #(t,t0)  propagates  the  state  from  time  to  to  time  t.  Another  valid 
expression  of  the  state  transition  matrix  is  the  differential  equation 

$(Mo)  =  A(t)$(t,t0)  (2.80) 

From  Reference  [15],  an  identity  of  the  state  transition  matrix  is 


$(*2,  to)  =  $(*2,  ti)#(ti,  t0)  (2.81) 

This  identity  is  useful  in  de-bugging  and  validation  of  code  for  the  estimator  program 
that  will  be  mentioned  in  the  next  chapter. 

2-4  Estimation  Theory 

Thus  far,  the  theory  governing  the  deterministic  dynamics  of  the  Target- 
Observer  system  has  been  discussed.  Following  the  detailed  derivation  provided 
by  Wiesel  [15],  this  section  will  briefly  explore  the  theory  behind  the  least  squares 
estimation  method  of  data  reduction. 

The  method  of  least  squares  estimation  is  based  primarily  on  the  following 
assumptions: 

•  Gaussian  distribution  of  error  in  the  observational  data 
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•  Deterministic  dynamics  (i.e.  no  error  in  the  dynamics  model) 

•  Principle  of  Maximum  Likelihood  approach  to  probability 

2-4-1  Probability  Theory.  French  mathematician  Abraham  de  Moivre 
(1667-1754)  first  introduced  the  mathematical  equation  of  the  normal  distribution 
curve  in  1733  [21].  The  normal  distribution  is  often  referred  to  as  the  Gaussian 
distribution  in  honor  of  the  German  mathematician  and  astronomer  Carl  Friedrich 
Gauss  (1777-1855).  Gauss,  considered  a  mathematical  Wunderkind  or  child  prodigy, 
is  credited  for  the  development  of  the  theoretical  basis  of  least  squares  in  1795  at  the 
age  of  eighteen  [22,  23].  Gauss  introduced  and  applied  his  method  of  least  squares  in 
1801  to  accurately  re-discover  the  asteroid  Ceres  (although  the  method  was  officially 
published  by  Adrien-Marie  Legendre  five  years  later)  [22,  23].  Beginning  with  the 
Gaussian  or  normal  statistical  distribution  model  [21],  explained  below,  the  familiar 
probability  density  function  is 

p(w) = w^exp{(^L)  (2'82) 

Following  the  derivation  provided  by  Wiesel  [15],  in  least  squares  estimation  theory, 
the  observational  data  (i.e.  the  raw  physical  measurements)  contain  error,  which  has 
a  Gaussian  distribution  expressed  as 

1  ( — e? 

where  P(e)  is  the  probability  of  obtaining  an  error,  e.  Once  a  range  of  error  is 
defined,  the  area  under  the  Gaussian  curve  will  give  the  probability  of  the  error 
occurring  within  the  range.  Lets  describe  a  measurement,  x,  as 

x  =  xq  +  e  +  b  or  e  =  x  —  b  —  x  o  (2.84) 


P(e)  = 


(2.83) 
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where  xq  is  the  unknown  true  value,  the  true  random  error  e  in  the  measurement, 
and  the  bias  b  of  the  instrument.  From  another  perspective,  applying  Equation  2.83 
and  the  error  function  to  observational  data,  the  probability  of  data  x  lying  within 
a  defined  interval  can  then  be  predicted.  The  distribution  is  normalized  so  that  the 
sum  over  all  values  gives  a  probability  of  1.  The  standard  deviation,  a ,  defines  the 
width  of  the  Gaussian  curve  which  is  symmetric  about  its  mean,  /i,  often  defined  as 
the  true  value  and  with  a  bias  of  zero  (no  shift).  A  large  a  results  in  a  broad  Gaussian 
curve  (i.e.  low  precision)  and  a  small  cr  results  in  a  narrow  curve  (i.e.  high  precision). 
Therefore,  an  instrument  is  considered  precise  when  it  has  a  relatively  small  standard 
deviation  because  the  measurement  data  points  have  a  smaller  spread  from  the  true 
value,  /.(.  The  Gaussian  curve  shown  in  Figure  2.5  illustrates  that  approximately 
68%  of  the  data  sample  (i.e.  measurements  x)  will  lie  within  the  range  of  ±1  a  from 
the  actual  true  value,  p,  =  x0. 


PM 


-3a  -2(7  -1(7  p  J  cr  2  a  3a 

Figure  2.5  Gaussian  Distribution,  (courtesy  of  Beirne) 


Before  proceeding,  the  Central  Limit  Theorem  should  be  mentioned.  Wiesel  [15] 
describes  this  theorem  as  “the  keystone  of  estimation  theory”  because  this  theorem 
is  what  allowed  Gauss  to  justify  the  assumption  that  the  statistical  behavior  of  error 
inherent  in  measurement  data  is  Gaussian  in  nature.  Although  the  lengthy  proof  is 
not  shown  here,  it  can  be  found  in  great  detail  in  the  text  by  Wiesel  [15]. 

In  the  simplest  case,  multiple  measurements  are  made  in  an  attempt  to  accu¬ 
rately  determine  a  single  true  value  that  one  wishes  to  know  exactly.  Unfortunately, 
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the  true  value  is  never  known  exactly  to  infinite  point  precision;  instead  one  must 
settle  for  an  estimate  of  the  true  value.  Given  N  independent  measurements,  aq,  of 
some  exact  quantity,  xo ,  the  joint  probability  is  given  by  the  product  of  the  individual 
Gaussian  distributions 


P(xi,x2,  ...,xN)  =  (27 r) 


AT 

2 


(Xj  -  £q)2\ 

2^2  ) 


(2.85) 


where  the  instrument  standard  deviation,  a,  describing  the  measurement  error  is 
known.  Keep  in  mind  the  accuracy  of  the  measurement  data  described  by  its  a- value 
is  actually  a  property  of  the  instrument  making  the  measurement.  As  mentioned  ear¬ 
lier,  since  the  true  value  that  one  attempts  to  know  through  multiple  measurements 
is  actually  unattainable,  and  the  true  error  of  the  measurement  is  also  unknown, 
Equation  2.85  leads  to  an  indeterminable  condition.  In  order  to  make  progress  in 
quantifying  an  unknown  quantity,  an  estimate  x  of  the  true  quantity  x.(,  will  be 
made.  The  pessimistic  and  ultra-conservative  approach  is  to  assume  that  the  error 
in  the  measurement  data  is  present  in  its  most  maximum  form  possible.  This  would 
lead  to  the  estimate  of  x  =  ±oo  which  would  lead  nowhere  in  estimating  the  true 
value.  Therefore,  the  optimistic  approach  or  Principle  of  Maximum  Likelihood  is 
used.  This  principle  defines  the  “...estimate  x  as  the  value  of  x0  which  maximizes 
the  probability  of  having  obtained  the  actual  data  set  [15].”  Equation  2.85  is  used 
with  the  substitution  of  the  estimate  x  for  the  true  value  x.q 


P  (xt)  =  (2tt)-t 


N 


) 


(2.86) 


The  next  step,  using  the  Principle  of  Maximum  Likelihood,  is  to  maximize  the 
probability.  By  taking  the  derivative  of  the  exponential  quantity  of  Equation  2.86 
with  respect  to  the  estimate  x  and  setting  it  equal  to  zero  minimizes  the  exponential 
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quantity  which  maximizes  the  probability  resulting  in  the  Method  of  Least  Squares 


d  ^ '  (xi  —  xf 


dx  ^  2  of 

1=1  1 


=  0 


(2.87) 


The  least  squares  method  is  commonly  described  as  minimizing  the  sum  of  the  square 
of  the  residuals.  Taking  the  derivative  yields 


=0 


2=1 


a ? 


(2.88) 


The  key  step  in  this  entire  formulation  was  replacing  the  unattainable  true  error  of 


e;  =  Xi  -  x0 


(2.89) 


with  the  residual 

Ti  =  Xi  —  x  (2.90) 

where  the  residual  is  dehned  as  the  difference  or  essentially  the  error  between  the 
observed  quantity  and  the  calculated  quantity.  Therefore,  the  statistics  of  the  resid¬ 
uals  can  give  you  insight  on  the  estimate.  If  the  estimate  is  nearly  equal  to  the  true 
value  then  the  residual  will  be  nearly  equal  to  the  true  error  [15]. 

2.4.2  Linear  Least  Squares.  Following  the  text  by  Wiesel  [15],  the  objective 
is  to  estimate  the  state  of  a  linear  dynamic  system  at  epoch  time  to  from  measurement 
data  z i  taken  at  specific  observation  times  tt . 


Xi(t)  =  &{ti,to)X{t0)  (2.91) 

In  linear  least  squares  the  key  assumption  is  that  the  system  state  at  time  tt  is 
linearly  related  to  the  data  observed  at  2*.  That  is,  the  data  taken  at  a  specific  time 
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z i{ti)  can  be  expressed  as  the  observation  relation 


z  i{ti)  =  HiX(ti)  +  e,( 


(2.92) 


where  the  data  is  simply  a  scalar  multiple,  Hi ,  of  the  state  X(f;)  and  the  true  error 
of  the  data  et.  By  substituting  the  dynamics  of  Equation  2.91  into  Equation  2.92 
relates  the  error  from  a  measurement  taken  at  a  current  time  to  the  state  at  epoch 
time  to 

=  Zi(ti)  -  Hi$(ti,t0)X(t0)  (2.93) 

where 

Ti  =  H^(tht0)  (2.94) 


Since  there  are  N  observations  it  is  convenient  to  consolidate  the  data  in  matrix 
format.  The  single  observational  data  vector  z,  may  actually  contain  more  than  one 
measurement  taken  at  a  single  time  tt  (e.g.  range,  azimuth,  elevation  data  taken 
simultaneously  at  each  time).  All  the  data  vectors  are  assembled  into  a  larger  data 
vector  spanning  t\  to 

Zl 
Zo 

>  (2.95) 

Z  N 

/ 

The  observation  matrix  is  formed  as 
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It  will  be  assumed  that  each  data  measurement  z,  is  statistically  independent  and 
therefore,  the  instrument  covariance  matrix  can  be  formed  of  individual  covariances 
on  the  diagonal  with  zeros  on  the  off-diagonals.  Each  covariance  Qi  corresponds  to 
each  observation  z, 


where  cq  is  the  standard  deviation  of  the  measurement  by  the  instrument,  also  de¬ 
noted  as  cTinstr  ■  The  instrument  covariance  matrix,  Q,  is  square  with  the  dimensions 
of  NxN  where  N  is  the  number  of  observations.  Resulting  from  the  derivation  pro¬ 
vided  by  Wiesel  [15]  the  estimate  of  the  linear  dynamical  system  X  at  epoch  time  to 
is  given  by 

X(t0)  =  (rTQ-1T)-1TTg-1z  (2.98) 

with  an  associated  state  covariance  Px  at  epoch  time  to 

Py(to)  =  (TtQ~1T)~1  (2.99) 

2-4-3  Nonlinear  Least  Squares.  With  the  linear  least  squares  method  de¬ 
fined,  the  nonlinear  method  will  now  be  developed  following  the  derivation  provided 
by  Wiesel  [15].  Nonlinear  least  squares  estimation  is  used  in  the  case  of  systems 
whose  dynamics  is  non-linearly  related  to  the  observations  made.  This  is  applicable 
to  most  real-world  systems  inherently  exhibiting  nonlinear  behavior.  Since  the  state 
matrix  of  the  system  in  this  thesis  includes  nonlinear  dynamics,  where  no  linear  rela¬ 
tionship  exists  between  the  states  and  the  observations,  this  nonlinear  least  squares 
method  will  be  used. 
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Begin  with  the  set  of  equations  of  motion  of  the  system  whose  solutions  will 


be  found  by  numerical  integration 


X  =  f(X,t) 

This  numerical  integration  would  result  in  the  explicit  solution  as  a  function  of  time 
and  the  initial  conditions  or  initial  states 

X(f)  =  h(X(f0)T)  (2.100) 

Developed  in  the  previous  section  of  Linearized  Dynamics,  the  linearization  of  a 
nonlinear  system  results  in  following  statement 

SX(t)  =  Q(t,t0)SX(t0) 

with  i)X  as  a  small  change  in  the  state,  and  the  linearized  dynamics  formed  into 
the  State  Transition  Matrix,  <h,  linearized  about  some  dynamical  trajectory.  The 
objective  is  still  to  know  the  true  state,  Xo,  which  is  unattainable.  Earlier  in  the 
subsection  of  Probability  Theory,  it  was  discussed  that  one  must  settle  for  an  estimate 
of  the  state,  X.  However,  to  initiate  the  nonlinear  least  squares  estimation  process 
of  finding  such  an  estimate,  one  must  initially  use  a  reference  trajectory  denoted  as 
Xre/.  Xre/  is  initially  considered  the  a  priori  estimate  of  the  initial  states  of  the 
dynamical  system  at  some  epoch  of  time  (i.e.  the  initial  educated  guess).  Using 
this  least  squares  method,  changes  in  the  state  <5X  will  be  generated  to  correct 
the  reference  trajectory  to  a  point  of  satisfaction  where  it  can  be  declared  that  the 
reference  trajectory  Xre/  is  the  estimated  trajectory  X  of  the  true  state  X0. 

The  observation  relation  as  a  function  of  the  state  variables  and  the  observation 
geometry  is  defined  as 


z  i{ti)  =  G  (X(ti),ti) 
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(2.101) 


The  true  error  in  the  observational  data  is  given  by  the  difference  between  the  actual 
collected  measurement  and  what  would  be  the  perfect  or  true  measurement. 


z  -  z0 

(2.102) 

G(X,  t)  —  G(X0,  t) 

(2.103) 

G(X0  +  SX,  t)  —  G(X0,  t) 

(2.104) 

dl6x w 

(2.105) 

As  mentioned  earlier,  it  is  expected  that  the  true  error  will  be  approximately  equal 
to  the  residual.  In  order  to  develop  a  relationship  describing  how  the  observation 
function  would  change  with  respect  to  a  change  in  the  system  states,  the  observation 
function  is  linearized  with  respect  to  the  state  variable  about  the  reference  trajectory 

(2.106) 

ref  (T  )  J’i 


*<*>-§ 


The  residual  vector  of  the  observation  can  be  calculated  by  the  difference  of  the  ac¬ 
tual  measurement  and  the  expected  measurement  defined  by  the  observation  Equa¬ 
tion  2.101 

r*  =  Zj  —  G(Xref(ti),ti)  (2.107) 

The  residuals  from  the  observation  data  with  Equation  2.94  follow  as 

r,  «  HiSX(U )  (2.108) 

=  (2.109) 

=  TiSX(to)  (2.110) 

Similar  to  the  linear  least  squares  method,  the  result  is 

sx(t0)  =  (rTg-1r)-1TTg-1r  (2.111) 
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Therefore,  the  equations  for  nonlinear  least  squares  estimation  are 


Psx(to)  =  {TT  Q~lT)~l  (2.112) 

where  after  convergence  would  be  defined  as 

Px  =  Psx  (2.113) 

and  the  estimated  trajectory  would  then  be  expressed  as 

X(t0)  =  Xre/(t0)  +  SX(t0)  (2.114) 

When  or  if  the  reference  trajectory  converges  to  a  solution  these  previous  equations 
describe  the  estimated  state  of  the  true  system  at  some  epoch  fo¬ 
il  5  Summary 

As  mentioned  earlier,  this  research  is  focused  on  estimating  the  relative  motion 
dynamics,  rotational  dynamics,  and  moments  of  inertia  of  the  RSO.  The  theory 
provided  in  this  chapter  yields  the  equations  of  motion  and  associated  assumptions 
governing  the  dynamics  of  the  Observer- Target  system,  estimation  theory,  and  the 
definition  of  the  state  vectors. 

For  relative  motion  dynamics,  the  equations  of  motion  for  relative  position 
(Equ.  2.23-2.25)  and  relative  velocity  (Equ.  2.27-2.29)  of  the  Observer- Target  system 
are  developed  and  consolidated  in  the  form  of  a  state  transition  matrix  (Equ.  2.31). 
The  equations  of  motion  for  the  moments  of  inertia  (Equ.  2.66-2.68)  are  provided 
as  well. 
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Six  ordinary  differential  equations  of  motion  are  required  to  describe  the  ro¬ 
tational  dynamics  of  the  RSO  rigid-body.  Euler’s  rotational  equations  of  motion 
(Equ.  2.43-2.45)  are  derived  along  with  the  remaining  three  complimentary  Euler 
rate  equations  (Equ.  2.60-2.62)  using  a  3-2-1  Euler  rotation  sequence. 

The  estimation  theory  discussed  in  this  chapter  covers  the  linearization  meth¬ 
ods  needed  for  the  dynamics  and  observation  function  and  the  mathematics  for 
implementing  the  non-linear  least  squares  method  in  a  computer  algorithm.  Also  of 
great  importance,  the  states  which  characterize  the  Observer-Target  system  at  any 
instance  in  time  is  provided  in  the  form  of  two  state  vectors  (Equ.  2.65  and  2.70). 
The  12-state  vector  is  used  for  estimation  cases  involving  relative  motion  and  rota¬ 
tional  dynamics  of  the  RSO.  In  estimating  the  rotational  dynamics  and  moments  of 
inertia  of  the  RSO,  the  9-state  vector  version  is  used. 

The  general  theory  and  assumptions  discussed  here  serve  as  the  basis  for  the 
approach  taken  in  the  effort  to  explore  the  estimation  of  the  dynamics  and  physical 
properties  (i.e  MOI)  of  the  RSO  using  range  observational  data.  The  theory  of  this 
chapter  is  implemented  in  the  methodology  developed  in  the  next  chapter,  specific 
to  this  research  effort. 


2-31 


III.  Methodology 

This  chapter  will  begin  by  formulating  the  observation  function  based  on  the  obser¬ 
vation  geometry  of  the  Observer- Target  system  involved  in  the  proximity  operation. 
The  formulation  of  the  linearized  observation  function  and  dynamics  will  follow.  A 
computer  program  written  in  MATLAB  was  created  to  model  the  two-body  dynam¬ 
ics  of  the  Observer-Target  system.  This  program,  termed  the  truth  model ,  is  used  as 
the  dynamics  model  of  the  system  as  well  as  a  generater  for  simulated  range  mea¬ 
surement  data.  Finally,  the  algorithm  for  the  nonlinear  least  squares  estimator  is 
discussed.  The  estimator  program,  written  in  MATLAB,  utilizes  the  simulated  range 
data  to  estimate  the  states  of  the  RSO  which  will  provide  insight  to  its  dynamical 
and  physical  properties. 

The  original  effort  of  this  thesis  was  to  explore  the  performance  of  a  nonlinear 
least  squares  filter  for  estimating  the  dynamics  of  a  Target.  The  rotational  velocities 
and  Euler  orientation  angles  were  of  particular  interest.  This  led  to  the  use  of  the 
12-state  dynamics  vector  as  developed  in  Equation  2.65.  Estimated  states  include 
the  relative  position,  relative  velocity,  angular  velocity,  and  Euler  angles.  In  an  effort 
to  explore  the  feasibility  of  estimating  the  moments  of  inertia  of  the  RSO  to  provide 
insight  into  the  object’s  physical  properties,  a  9-state  version  of  the  state  vector 
was  created  as  shown  in  Equation  2.70.  In  a  scenario  where  the  relative  position 
and  velocities  are  considered  known  to  sufficient  accuracy,  the  estimated  parameters 
include  the  angular  velocity,  Euler  angles,  and  moment-of-inertia  components  of  the 
RSO.  This  chapter  will  focus  on  the  development  of  the  H  and  $  matrices  pertaining 
to  the  12-state  variable  version.  The  formulation  method  for  the  9- variable  version 
is  identical  and  the  program  algorithms  are  not  affected.  Details  specific  to  the 
9- vector  version  are  found  in  the  Appendix  A. 
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3.1  Observation  Function 


During  the  space  proximity  operation,  the  Observer  spacecraft  will  be  collect¬ 
ing  range  measurements  of  continuously  tracked  points  on  the  Target  body  for  the 
duration  of  the  observational  period.  The  relative  position  and  velocity  states  of 
the  Observer  as  well  as  the  Target’s  rotational  states  at  a  given  instant  in  time  are 
related  to  the  range  observation  made  at  that  particular  time.  This  leads  to  the 
purpose  of  the  observation  function.  The  observation  function  is  an  expression  that 
is  based  on  the  observation  geometry  shown  in  Figure  3.1  that  results  in  expected 
observations  as  a  function  of  the  dynamic  state  variables.  Therefore,  given  the  dy¬ 
namic  states  of  the  Observer- Target  system,  the  range  measurement  that  should  be 
made  by  the  Observer  spacecraft  can  be  predicted.  This  predicted  range  measure¬ 
ment  and  the  actual  range  measurement  made  (using  simulated  range  data)  is  at 
the  heart  of  the  least-squares  estimation  process.  The  observation  relation  is  derived 


&N 


Figure  3.1  Observation  Geometry  of  Observer- Target  System 
starting  with  the  relative  position  vector  8p.  Careful  attention  is  paid  to  the  starting 
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and  terminating  points  of  the  vectors. 


5  f)  Yn 


(3.1) 


solving  for  rn,  yields 


rn  =  Rn  -  8p 


(3.2) 


where  the  subscript  n  denotes  observation  point  n  on  the  Target  body.  The  vector 
Rn  represents  the  position  of  point  n  from  the  center-of-mass  of  the  Target  in 
The  vector  rn  is  representative  of  the  observation  vector  as  seen  from  the  Observer 
spacecraft  to  the  observational  point  on  the  Target  body  in  $e.  Although  the  rn 
vector,  by  definition,  contains  magnitude  (i.e.  range)  and  directional  information, 
this  thesis  will  only  explore  range  observational  data.  Therefore,  the  actual  obser¬ 
vational  quantity  is  simply  the  magnitude  of  the  rn  vector,  yielding  range ,  denoted 
by  the  scalar  rn 

rn  =  ||r„||  =  ||Rn  -  Sp\\  (3.3) 

However,  note  that  the  observation  function  contains  vectors  of  mixed  basis.  Per¬ 
forming  a  coordinate  frame  transformation  to  $e  using  the  proper  rotation  matrix 
(Equ.  2.54)  with  the  previous  equation  results  in  the  expression  below  with  its  time- 
dependence  shown 


r:  = 


kwh  = 


.[R^R'-ip'll 

K  -  Sp‘\\ 

ir;w- w)ii 


(3.4) 

(3.5) 

(3.6) 
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define  the  position  vectors  of  the  nth  observation  point  as 


M 

Pxn(t)  j 

b%,  >  and  RJ(t)  =  < 

Pyjt)  t 

UJ 

Vznit)  J 

(3.7) 


Substituting  Equation  3.7  for  R®  and  Equation  2.4  into  Equation  3.6  yields 


rnll  =  \\(PxneR+pynee+pZneN)  -  (5xeR  +  5yee  +  5zeN)\\  (3.8) 

rnll  =  II  fen  -  Sx)eR  +  fe„  -  5y)ee  +  fen  -  ^)ejv)||  (3.9) 

rn  =  \! fen  -  fe)2  +  fen  -  fry)2  +  fen  “  (3-10) 


Finally,  the  range  observation  function  G  is  defined  below  in  terms  of  the  dynamic 
state  variables  X  which  will  predict  the  range  measurement  data.  This  function,  as 
developed  here,  is  unique  to  the  system  geometry  of  this  research  as  illustrated  in 
Figure  3.1. 


G  =  range  =  yj fen  -  fe)2  +  fen  -  ~5y)2  +  fen  -  Jz)2  (3.11) 

In  the  case  of  multiple  range  measurements  to  multiple  points  taken  essentially 
simultaneously,  the  function  G  simply  becomes  a  vector  of  range  measurements 
for  a  single  observation  time 


/  \ 

ri 

V fei  -  +  fei  -  <fe2  +  fei  -  fe)2 

G  =  < 

r2 

>  =  < 

V fe2  -  Sx)2  +  fe2  -  Sy )2  +  fe2  -  fe)2  ^ 

rN 

V  / 

k  Vfeiv  - Sx)2  +  (pvn  -  sv)2  +  fejv  -  fe>2  , 

for  observation  points  1  to  AT;  that  is,  n  =  1,2,...  ,7V. 
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3.2  Linearization  of  Observation 

This  section  is  devoted  to  linearizing  the  observation  function  to  find  a  linear 
relationship  between  the  change  of  the  observation  with  a  change  in  the  dynamics 
state  vector.  This  expression  will  be  denoted  as  the  H  matrix  and  will  be  needed 
for  the  estimation  algorithm. 

In  terms  of  the  state  variables,  note  that  the  final  version  of  the  range  observa¬ 
tion  expression,  Equation  3.11,  is  a  function  of  the  CW  relative  position  vector  and 
the  position  vector  of  the  observation  point  in  $e.  However,  the  observation  point 
in  is  a  function  of  the  Euler  orientation  angles  due  to  the  rotation  of  the  body 
frame  of  the  Target.  The  version  of  the  observation  function  in  terms  of  the  relative 
position  state  variable,  bp,  and  the  Target  observation  point  position,  a  non-state 
variable,  is  given  by  Equation  3.11  as 

G  =  rn  =  y l\ (pXn  -  Sx )2  +  ( pVn  -  by)2  +  (p*n  -  Sz )2 

From  Equation  2.106,  recall  the  linearized  observation  function  with  respect  to  the 
state  vector  is  defined  as 


Ui(ti)  = —(Xref(ti),ti)  (3.13) 

Note  that  the  range  equation,  Equation  3.11,  is  a  function  of  Sx,  by,  bz,  pXn ,  pVn , 
and  pZn  at  any  moment  in  time.  In  terms  of  the  dynamics  state  variables,  the  range 
equation  is  only  a  function  of  Sx,  by,  bz,  0\ ,  02,  and  63.  For  simplicity,  all  other  state 
variables  which  the  range  equation  is  not  a  function  of  are  ignored  till  the  end.  This 
results  in 


H  = 


dG 

dX 


dG  dG  dG  dG  dG  dG 

8Sx  dSy  85z  86\  862  863 


(3.14) 
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Recalling  that  pXn ,  pVn ,  and  pZn  are  functions  of  the  Euler  angles  at  specific  moments 
in  time  due  to  the  rotation  of  the  Target  body,  the  H  matrix  can  simply  be  found 
using  the  chain  rule  in  the  following  fashion 


(315) 


The  partial  derivatives  shown  above  (and  found  in  Appendix  A)  are  coded  into  the 
program  and  the  matrix  multiplication  is  carried  out  by  the  computer  program.  The 
state  variables  which  are  not  involved  are  simply  zero.  This  results  in  the  solution 
of  H,  for  every  observation  time  tt. 

In  terms  of  the  12- state  vector,  defined  as 

r  -|  t 

X  =  Sx  Sy  5z  Sx  Sy  Sz  cl>2  w 3  62  #3 

the  linearized  observation  function  is  the  solution  to 


]J  _  dG  dG  dG 

d  Sx  dSy  dSz 


1 

0 

1  0 

dG  dG  dG 

dpxn  dpVn  dpZn  \  q 

0 

0 


0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

dPxn 

dPx„ 

dpx„ 

ddi 

d02 

de3 

0 

0 

dpyn 

dpyn 

dpyn 

d6i 

de2 

d6  3 

0 

0 

dpzn 

dpzn 

9Pzn 

<96>1 

de2 

d03  . 

H 


qg  8G  8G  o  n  n  n  n  n  Qg  og  qg 

dSx  d8y  dSz  U  U  U  U  U  U  d9 1  d62  dd3 


(3.16) 


In  terms  of  the  9- state  vector,  defined  as 


X  —  Co-'x  ui2  ^3  9\  62  O3  A  B  C 


where  the  relative  position  variables  are  considered  known  constants  for  each  obser¬ 
vation  time  and  are  then  not  estimated  states,  the  linearized  observation  relation  is 


3-6 


simply  the  solution  to 


H 


0  0  0  || 


dG 

oe  2 


dG 

ao3 


0  0  0 


(3.17) 


3. 3  Linearization  of  Dynamics 

This  section  will  show  how  the  dynamics  of  the  system  is  linearized  resulting  in 
the  formulation  of  the  state  transition  matrix,  <E.  This  formulation  will  be  shown  for 
the  12-state  vector  version.  The  9-state  vector  version  will  simply  follow  the  same 
method  and  is  given  in  Appendix  B. 

Recall  from  Equation  2.31  that  the  CW  state  transition  matrix  dy,,,  is  already 
known.  Therefore,  the  state  transition  matrix  for  the  rotational  portion  of  the  dy¬ 
namics  is  in  need.  As  defined  in  Chapter  2,  the  rotational  state  vector  is 


OJl  U>2  Ll>3  6\  62  $3 


Therefore,  applying  Equation  2.71,  the  time-derivative  of  the  state  vector  is  shown 
below.  The  vector  function,  f ,  is  then  associated  with  each  state  equation  of  motion. 
Functions  /1,  /2,  /3,  correspond  to  Equations  2.43,  2.44,  2.45  and  functions  /4 ,  /5, 
/6,  correspond  to  Equations  2.60,  2.61,  and  2.62,  respectively. 


Xrot(t)  — 


\ 

<hi 

/  \ 

h 

Ll>2 

/2 

h  3 

>  =  < 

h 

0i 

U 

$2 

h 

e3 

\  / 

=  m 


(3.18) 
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From  Equation  2.77,  the  partial  derivatives  of  the  equations  of  motion  with  respect 
to  the  state  variables  of  the  rotational  dynamics  are  symbolically  shown  below  (the 
solution  is  found  in  Appendix  B) 

dfi  dfi  dfi  dfi  dfi  dfi 
dun  8u2  du>3  801  80 2  803 

d/2  d/2  d/2  d/2  d/2  d/2 
duo  1  duo 2  duo 3  d#i  d#2  d03 

d/3  d/3  d/3  dfs  dfs  dfs 
duo  1  da;2  du;3  d#i  d#2  d03 
d/4  d/4  d/4  d/4  d/4  d/4 

da;i  do;2  du;3  d6>i  d<92  d<93 

d/5  d/s  d/5  d/5  d/5  d/5 

do;i  da;2  du;3  d#i  d#2  d#3 
d/e  d/6  d/6  d/6  d/6  d/6 
do;i  do;2  du;3  d#i  d#2  d03 

The  state  transition  matrix  for  the  rotational  portion  is  defined  as 


011  012  013  014  015  016 

021  022  023  024  025  026 

T  031  032  033  034  035  036 

$rot  = 

041  042  043  044  045  046 

051  052  053  054  055  056 

061  062  063  064  065  066 

and  the  time-derivative  of  the  state  transition  matrix  is  simply  defined  as 


(3.20) 


(3.19) 


011 

012 

013 

014 

015 

016 

021 

031 

022 

032 

023 

033 

024 

034 

025 

035 

026 

036 

041 

042 

043 

044 

045 

046 

051 

052 

053 

054 

055 

056 

061 

062 

063 

064 

065 

066 

(3.21) 
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The  matrices  above  are  coded  into  the  estimation  program  and  using  a  Runga-Kutta 
algorithm,  the  numerical  solution  to  Equation  2.80 


^rot(G^o)  A(i)i&rot (G  toi) 


is  solved  for  resulting  in  &rot(t,  t0 ).  Incorporating  the  CW  portion  (Equ.2.31),  results 
in  the  newly  developed  12  x  12  state  transition  matrix  for  the  12-state  vector  version 
used  in  the  estimation  process 


$(Mo) 


a>6x6 

CW 

q6x6 

q6x6 

<i>6x6 

rot 

(3.22) 


Following  the  same  methodology,  with  the  omission  of  the  CW  states  and  changes 
to  the  A  matrix  to  include  the  moment-of- inertia  terms,  the  9-state  vector  version 
of  the  9x9  state  transition  matrix  is  developed  and  shown  in  Appendix  B. 


3.4  Truth  Model  and  Data  Generator 

A  truth  model  of  the  dynamics  was  created  to  provide  a  basis  for  comparison 
between  the  simulated  actual  dynamics  and  the  estimated  dynamics.  Also,  the  truth 
model  provides  a  method  of  generating  observational  range  data  since  no  real-world 
experimental  data  was  available  for  this  effort.  Both  the  Truth  Model  and  the  Data 
Generator  were  coded  in  MATLAB. 

The  Truth  Model  generates  two-body  deterministic  dynamics  data  for  relative 
position,  relative  velocity,  Target  angular  velocity,  Target  Euler  angle  orientation, 
and  Target  moments  of  inertia  for  the  Observer-Target  system  for  the  duration  of 
the  proximity  operation  in  which  range  observations  are  being  made.  Assumptions 
programmed  into  the  Truth  Model  are  as  follows: 

•  Target  RSO  is  in  a  true  circular  orbit  about  the  Earth 

•  Moments  of  Inertia  of  the  Target  are  constant 
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•  No  perturbations  exist 

The  Data  Generator  is  a  program  which  uses  the  Truth  Model  to  provide  a  dy¬ 
namical  scenario  in  which  observational  range  data  would  be  taken  by  the  Observer 
spacecraft.  With  the  dynamics  data  generated  for  each  moment  in  time,  the  obser¬ 
vation  function  (Equ.  3.11)  is  used  to  generate  observational  range  measurements. 
Two  forms  of  range  measurements  are  created  by  this  program. 

•  Perfect  range  measurements  with  no  error 

•  Noisy  range  measurements  containing  Gaussian  noise 

The  pseudo-random  Gaussian  noise  is  added  to  the  perfect  range  data  to  simulate 
real-world  data  containing  measurement  error.  Assumptions  incorporated  into  the 
Data  Generator  are  as  follows: 

•  The  time  vector  associated  with  the  dynamics  propagation  is  identical  to  the 
observation  time  vector. 

•  Observation  measurements  are  taken  at  constant  time  steps.  A  non-continuous 
data  scenario  is  not  considered. 

•  Observational  range  sensor  and  pattern-recognition  capability  provides  an  abil¬ 
ity  to  track  user-defined  points  on  the  Target  body  for  the  duration  of  the  ob¬ 
servation  time  vector.  Data  acquisition  limitations  due  to  physical  obstructions 
are  not  considered. 

•  For  the  simulated  range  data  with  noise,  the  standard  deviation  of  instrument 
error  making  the  range  measurements  is  known  and  is  constant  through  the 
observational  period. 

•  All  measurements  are  statistically  independent 

•  Target  observation  points  are  known  in 

The  program  architecture  of  the  Truth  Model  (including  the  Data  Generator) 

developed  for  this  thesis  is  shown  in  Figure  3.2. 
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Truth  Model  Program  Architecture 


CW  Equations  of  Motion 


-  Evaluate  CW  position  solutions 

-  Evaluate  CW  velocity  solutions 

-  Xe4t) _ 


Numerical  Integ  of  Angular  Velocity  ODEs 
Numerical  Integ  of  Euler  Rate  ODEs 
■  X.ptft) _ 


ECI  Orbit  Calculator 


Orbit  of  Target  in  ECI-frame 
Orbit  of  Observ  in  ECI-frame 


3-1-3  Euler  Rotation  [R’e 
e-frame  to  ECI-frame 


3-2-1  Euler  Rotation  [Reb] 
b-frame  to  e-frame 


Generate  White  Gaussian  Noise 
Generates  Range  Data  w/o  Gaussian  Noise 
Generates  Range  Data  w/  Gaussian  Noise 


USER  INPUT  SCRIPT 


-  Time  Range  for  Dyn  &  Observations 

-  Dynamics  12-State  Vector  at  Epoch,  Xdyn(to) 

-  Moments  of  Inertia 

-  Orbital  Parameters 

-  Target  Observation  Point  Location 

-  Output  Settings 


TRUTH 


MODEL 


DYNAMICS 


-  Propagate  CW  Equations  of  Motion 

-  Progagate  Rotational  Equations  of  Motion 

-  Form  Dynamics  State  Model 


OBSERVATION 

-  Perform  Coordinate  Transformation 

-  Evaluate  Observ  Relation,  z  =  G(X) 

-  Form  Simulated  Range  Observ  Data 


DYNAMICS  STATE  DATA 


-  Simulated  12-State  Dyn  for  Time  Range 

-  Xdyn(t) 

-  Written  to  file 


RANGE  OBSERVATION  DATA 


-  Observation  Time  Vector 

-  Simulated  Range  Data  Vector  w/  Noise 

-  Simulated  Range  Data  Vector  w/o  Noise 

-  Written  to  three  seperate  files 


-  Various  Plots  of  Dynamics  &  Observations 

-  Gaussian  Noise  Histogram 

-  Observ  Geometry  Animation 

-  Print  Table  Data  to  Screen 


Target  Observation  Point  File 

-  Rn  =  \pXn  Pyn  Pz„\ 

-  Entered  by  the  User  in  matrix  format 

-  to  x  3  where  to  =number  of  points 

Figure  3.2  Truth  Model  Architecture 

3-4-1  Program  Execution.  Following  Figure  3.2,  in  the  user-input  script, 
the  user  enters  the  initial  values  for  the  state  variables  at  epoch.  That  is,  the  initial 
conditions  for  the  CW  relative  position,  CW  relative  velocity,  angular  velocity  of 
the  Target,  and  Euler  orientation  angles  of  the  Target,  as  well  as  moments  of  inertia 
are  defined.  The  altitude  of  the  Target  RSO  is  also  entered  as  well  as  the  epoch 
start  time,  time-step,  and  propagation  time  in  terms  of  orbital  period.  For  the 
observation  portion  of  the  Truth  Model,  the  number  of  Target  observation  points  to 
generate  simulated  range  data  for  is  entered.  The  standard  deviation  of  the  range 
data  is  also  entered.  A  separate  file  containing  the  vector  positions  of  the  observation 
points  in  3/  is  created  by  the  user.  This  provides  the  initial  conditions  and  constants 
needed  for  the  dynamics  portion  of  the  Truth  Model.  Using  the  CW  solutions  for 
relative  position  and  velocity,  the  CW  states  are  solved  for  every  time,  tt.  Using 
the  built-in  MATLAB  Runga-Kutta  ode45  routine,  the  numerical  solution  to  the 
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rotational  equations  of  motion  are  determined  for  tL .  This  results  in  the  propagated 
dynamics  which  is  written  to  a  file.  The  observation  portion  of  the  Truth  Model 
reads  the  dynamics  data  file  and  the  file  containing  the  observation  points,  along 
with  the  user’s  input  of  error  to  generate  simulated  range  observation  data  with 
noise  and  without.  The  range  data  along  with  the  observation  time  vector  is  written 
to  separate  files  for  use  by  the  Estimator  program.  Various  plots,  histograms,  and 
tables  are  also  produced. 

3-4-2  Program  Validation.  Several  methods  were  used  in  the  validation 
of  the  Truth  Model  code.  For  the  CW  relative  position  and  velocity  dynamics,  the 
analytical  solution  was  simply  verified  by  hand-calculations  for  arbitrary  data  points. 
Also,  test  cases  using  relative  velocity  components  in  the  cross-track,  radial-track, 
and  in-track  directions  were  used  to  check  for  the  expected  effects  of  orbital  drift, 
eccentricity  changes,  inclination  changes,  and  combinations  of  such. 

The  validation  of  the  rotational  dynamics  portion  of  the  model  was  more  in¬ 
volved.  Hand-calculations  were  performed  for  only  single-axis  rotations  to  verify 
angular  velocities  and  Euler  orientation  angles.  Also,  using  the  Target’s  moments  of 
inertia  and  angular  velocities,  test  cases  were  accomplished  to  check  for  rotational 
stability  about  the  major,  minor,  and  intermediate  axes.  Finally,  conservation  of 
angular  momentum  and  conservation  of  energy  checks  were  built  into  the  code  using 

||H||  =  ||lu;6i||  =  constant 

KE  =  -ojbl  ■  Iu>bt  =  constant 
2 

where  the  magnitude  of  the  angular  momentum  H  is  constant  in  (Equ.  2.34). 
The  kinetic  energy  of  the  rigid  Target  body  is  denoted  as  KE  and  given  by  the 
dot  product  shown  above.  These  checks  scan  through  the  generated  data  to  ensure 
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constant  values  of  angular  momentum  and  kinetic  energy  are  within  a  tolerance  set 

to  nr8. 

3.4.3  Simulated  Range  Data.  Using  the  vector  location  of  the  observation 
point  and  the  CW  relative  position  vector  from  the  generated  dynamics  data,  the 
observation  function  (Equ.  3.11)  is  evaluated  to  produce  simulated  range  data. 
Based  on  the  probability  theory  discussed  in  Chapter  2,  the  spacecraft  observation 
sensor  will  be  plagued  with  measurement  errors  following  a  Gaussian  distribution 
(Fig.  2.5);  Gaussian  noise  is  added  to  the  perfect  range  data  to  simulate  real-world 
data  containing  measurement  error.  A  MATLAB  pseudo-random  number  generator 
is  used  to  create  the  Gaussian  noise  which  is  initialized  with  every  execution  using  the 
current  computer  clock  time  as  the  seed.  Figure  3.3  illustrates  a  sample  distribution 
of  the  Gaussian  noise  for  an  arbitrary  data  set  spanning  one  orbital  period. 


Gaussian  Noise  Histogram 


Range  Error  Bin  Size  [m] 


Figure  3.3  Sample  Gaussian  Noise  Distribution 

Research  indicated  that  commercially  available  LIDAR  sensors  claim  range 
accuracies  from  1  m  to  7mm  (one  sigma)  from  distances  of  near  100  m  [24,  12,  25]. 
Since  this  thesis  is  focused  on  estimating  the  dynamical  properties  of  an  RSO  for  the 
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intention  of  performing  a  proximity  operation  possibly  involving  spacecraft-to-target 
interaction,  a  range  accuracy  under  one  meter  is  a  realistic  expectation.  Therefore, 


for  the  purpose  of  this  research,  the  range  measurement  accuracy  will  be  set  to  25  cm 
(one  sigma)  for  all  cases.  The  noise  error  simulated  in  a  sample  data  batch  is  shown 
in  Figure  3.4 

Range  Measurement  Error  from  Gaussian  Noise 
a  =  0.25  |r  =  0 
Target  Observation  Point  1 


-0.8 


4.1 - ' - ' - ' - ' - ' - 1 

0  1000  2000  3000  4000  5000  6000 

Time  [s] 


Figure  3.4  Sample  Range  Measurement  Noise  Error 


3.5  Nonlinear  Least  Squares  Estimator 

This  section  will  discuss  the  nonlinear  least  squares  estimator  that  was  created 
to  estimate  the  states  of  the  RSO.  This  program  written  in  MATLAB  is  illustrated 
in  Figure  3.5.  Two  different  approaches  of  the  estimation  are  attempted.  Although, 
both  versions  involve  the  same  process  and  same  methodology,  the  outcome  is  dif¬ 
ferent. 

The  first  estimation  approach  uses  the  standard  12-state  vector  (Equ.  2.65)  to 
estimate  the  relative  motion  and  rotational  dynamics  of  the  Target.  This  method 
provides  an  estimate  of  the  CW  position,  CW  velocities,  angular  velocities,  and 

Euler  orientation  angles.  The  second  approach  used  in  this  thesis  was  motivated 
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by  an  effort  to  explore  the  feasibility  of  estimating  the  moments  of  inertia  of  the 
RSO  directly  based  on  range  observations.  In  this  effort,  the  relative  position  data 
is  considered  known-processed  data  through  some  prior  estimation  filter  or  data 
refinement  method.  That  is,  the  relative  position  data  will  not  be  included  in  the 
least  squares  estimation  process;  rather,  it  will  be  supplied  to  the  algorithm  but 
treated  as  known  information.  In  this  9-state  version  (Equ.  2.70),  the  state  variables 
estimated  are  the  angular  velocities,  Euler  angles,  and  moments  of  inertia  of  the 
RSO.  This  method  allows  the  estimator  to  focus  its  effort  solely  on  estimating  the 
rotational  dynamics  and  moments  of  inertia  of  the  Target,  providing  insight  into 
both  its  dynamical  and  physical  properties.  The  results  of  both  approaches  will  be 
presented  in  Chapter  4. 


Estimator  Program  Architecture 


-  Form  complete  <F(ii,fo) 

-  12  or  9  state  vector  version 


-  Computes  &Cw{ti,to) 


-  Computes  $rot{U,to) 

-  incl  MOI  version 


-  Computes  Hi 

-  12  or  9  state  vector  version 


USER  INPUT  SCRIPT 


-  Estimated  Dynamics  State  Vector  at  Epoch 

-  Estimated  Moments  of  Inertia 

-  Instrument  a 

-  Filter  sensitivity  settings 

-  Data  output  settings 


N-LSQ  ESTIMATOR 


-  Propagate  CW  Equations  of  Motion 

-  Progagate  Rotational  Equations  of  Motion 

-  Form  Xref 

-  Compute  <F(q,fo) 

-  Compute  r,  =  z,  —  G»(X) 

-  Compute  Ti 

-  Compute  Running  Sums 

-  Calculate  Pgx 

-  Calculate  5X(to) 

-  Correct  Xref+1  =  Xre/(f0)  +  <SX(f0) 

-  Check  convergence  criteria 

-  Repeat  algorithm  as  necessary 


ESTIMATED  STATE  DATA  FILE 


-  X(i0) 

-X(t) 

-  Written  to  file 


COVARIANCE  OF  ESTIMATE 


Px 

-  Written  to  file 


-  Various  Plots  of  Estimated  Dynamics 

-  Histogram  of  Residuals 

-  Data  Range  and  Estimated  Range 

-  Print  Table  Data  to  Screen 


-  Rn  =  \pxn  Pyn  PzJ 

-  Entered  by  the  User  in  matrix  format 

-  m  x  3  where  m  =number  of  points 


-  Observation  time  vector 

-  Range  data  without  noise,  z 

-  Range  data  with  noise,  z 

Figure  3.5  Estimator  Architecture 
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3.5.1  Algorithm.  Following  the  theory  provided  in  Chapter  2  for  nonlinear 
least  squares  and  following  the  algorithm  by  Wiesel  [15],  the  estimation  algorithm  is 
as  follows: 


1.  Input  a  priori  estimate  of  states  at  epoch 

•  X(t0) 

2.  Read  data  files  produced  by  Truth  Model 

•  Observation  time  vector 

•  Range  data  with  or  without  noise 

•  File  containing  target  observation  points 

3.  Process  all  observational  data  for  each  time  t, 

•  Propagate  state  vector,  X(tj) 

•  Calculate  t0) 

•  Calculate  r*  =  z,  —  G(X) 

•  Calculate  Tj  =  0) 

•  Calculate  and  store  the  running  sum  'YhiT'[ Q^xTi 

•  Calculate  and  store  the  running  sum  JT  T^Q^Yi 

4.  Calculate  the  correction  covariance 

•  p^x  =  Qpr\ 

5.  Calculate  epoch  state  correction 

.  8X(t0)  =  hPsx^Ei^Q-1^ 

6.  Correct  the  Reference  Trajectory  at  epoch 

•  Xrep|_i(to)  =  Xref(to )  +  <5X(to) 
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7.  Determine  Convergence 


•  Covariance  criteria  check:  |£Xj(i0)|  TC  k2 \J Psxu 

•  Residuals  a  criteria  check:  crm.,tr  «  crresid,  tolerance  defined  by  k 3 

•  If  criteria  is  not  satisfied,  begin  next  iteration  at  step  3.  Otherwise, 
proceed  to  step  8 

8.  Declare  Estimate 

•  Estimate  trajectory  is  X(t0)  =  Xre/(t0) 

•  Covariance  of  estimate  is  Py  =  Psx 


3.5.2  Validation.  Several  methods  were  used  to  verify  the  correctness  and 
functionality  of  the  estimator  and  its  major  modules.  For  state  propagation  and 
observation  relation  calculations,  the  estimator  program  calls  the  same  validated 
script  files  used  by  the  validated  Truth  Model  program.  The  linearized  observation 
relation  was  check  by  this  approximation 


H, 


Gj(Xj  +  8Xj)  —  Gj(Xj) 


*X; 


(3.23) 


and  the  linearized  dynamics  was  checked  similarly  with  this  relation 


Hi.j 


f*(Xj  +  hXj)  —  fi(Xj) 


hX, 


(3.24) 


This  was  accomplished  for  rows  i  and  columns  (i.e.  state  variables)  j  using  small 
(~  10-4)  state  changes  3X?,  one  state  variable  at  a  time,  for  arbitrary  test  cases. 
The  state  transition  matrix  was  checked  to  ensure  this  property  was  upheld 


&(t2,to)  =  &(t2,ti)&(ti,t  0) 


(3.25) 
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Overall  functionality  of  the  estimator  was  checked  by  setting  the  a  priori  state 
estimates  at  t0  to  the  actual  true  states  at  t0  used  by  the  truth  model  to  generate 
the  data.  As  expected,  a  perfect  guess  generated  a  perfect  estimate  (i.e.  the  true 
trajectory)  with  zero  residuals  for  all  observations. 

3.5.3  RSS  Guess  Generator.  Step  1  of  the  estimation  algorithm  calls  for 
the  user  to  enter  the  a  priori  estimate  of  the  states  at  epoch,  X(f0)-  Simply  put,  the 
key  to  success  in  this  estimation  process  is  having  a  good  initial  guess  of  X(t0)  for 
the  estimator  to  begin  its  iterative  routine  with.  In  general,  there  is  no  guarantee 
that  the  algorithm  will  converge  on  the  actual  trajectory,  unless  the  initial  guess 
matches  the  actual  true  states  at  epoch  exactly.  But  an  initial  estimate  sufficiently 
close  to  the  truth  will  drastically  improve  the  outcome.  The  least  squares  algorithm, 
through  its  iterative  process,  strives  to  minimize  the  sum  of  the  residuals  squared, 
or  the  function  S 

S  =  rTQ-1r  (3.26) 

This  function  is  minimized  to  a  local  minimum.  However,  this  function  may  ac¬ 
tually  have  other  extrema  besides  the  local  minimum  where  the  desired  solution 
lies  [15].  Movement  to  any  other  local  minimum  would  result  in  a  case  where  the 
algorithm  does  not  converge  to  the  desired  solution.  For  this  reason,  the  Residual 
Sum  Squared  (RSS)  guess  program  was  created.  This  program  attempts  to  find  the 
global  minimum  using 

RSS(XuX2)  =  5>  -  G)2  (3.27) 

i 

as  a  function  of  two  state  variables  that  drive  S  — »  0.  It  should  be  noted  that  use 
of  the  RSS  guess  program  is  optional  and  independent  of  the  estimation  program. 

This  RSS  Guess  Generator  will  only  provide  a  best-guess  solution  for  two  state 

variables  of  interest.  It  is  assumed  that  the  user  has  a  priori  knowledge  of  all  other 

variables  with  a  certain  level  of  confidence  for  input  into  the  estimation  algorithm. 

The  limitation  to  two  state  variables  for  this  program  is  simply  due  to  the  multi- 
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dimensional  nature  of  this  problem  and  the  computational  intensity  involved.  In 
this  thesis,  the  rotational  states  of  the  Target  are  of  particular  interest.  Therefore, 
the  guess  generator  will  be  used  to  estimate  two  of  the  three  rotational  velocity 
components  of  the  Target.  The  RSS  guess  generator  algorithm  flows  as  such 

1.  Select  two  state  parameters  of  interest  for  a  generated  guess 

•  Choose  two  of  the  three:  u>i,  lo2,  uj. 3 

2.  Input  the  guess  bound  for  each  state  variable 

•  Parameter  1:  a  <  ui  <b 

•  Parameter  2:  c  <  ui  <  d 

3.  Input  the  guess  increment  or  step- size  within  the  bound 

•  Create  a  vector  of  parameter  1  values  from  a  to  b 

•  Create  a  vector  of  parameter  2  values  from  c  to  d 

4.  Compute  the  predicted  range  measurements 

•  Using  every  combination  of  parameter  1  and  parameter  2  values,  aug¬ 
mented  with  the  remainder  of  the  states,  assemble  X(f0),  and  propagate 
the  dynamics  to  each  observation  time  t,. 

•  Compute  G(X)  for  each  observation  time  tt 

•  Read  range  data  from  Truth  Model  output  file,  z, 

5.  Compute  RSS  solution 

•  Compute  sum  RSS  =  —  G)2  for  the  data  batch  as  a  function  of 

parameters  1  and  2 

6.  Output  best-guess  solution 
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•  Output:  3-D  surface  plot  where  the  rc-axis  corresponds  to  parameter  1  val¬ 
ues,  y-aixis  corresponds  to  parameter  2  values,  and  the  2- axis  corresponds 
to  the  RSS  value  for  the  combination. 

•  Locate  the  global  minimum  value  of  RSS. 

•  Using  the  RSS  global  minimum,  provide  the  corresponding  parameters  1 
and  2  values  as  the  best  guess  solution. 

If  computational  time  is  not  of  concern,  and  the  values  of  two  state  variables  at  epoch 
are  not  well  known,  this  tool  can  provide  an  extremely  good  estimate  for  the  two 
parameters  given  a  range  that  may  span  three  orders  of  magnitude.  Implementation 
of  this  tool  will  be  shown  in  the  next  chapter. 

3.5.4  Program  Execution.  To  begin,  the  user  enters  the  a  priori  estimate 
of  the  states  at  epoch  with  or  without  the  help  of  the  RSS  Guess  Generator.  It  is 
also  here  that  the  user  must  select  the  12-state  or  9- state  version  for  the  estimation 
process  and  the  number  of  observational  points  to  use.  Also  selected  is  the  data 
mode  for  the  processing  of  range  data  with  noise  or  without  noise.  The  user  has 
control  of  tuning  the  sensitivity  of  the  estimator  algorithm  by  adjusting  the  filter 
setting  for  the  k\,  k2 ,  and  k3  scaling  coefficients.  This  will  be  discussed  in  the 
next  subsection.  The  maximum  number  of  iterations  and  the  instrument  covariance 
value  in  terms  of  ainstr  are  also  defined.  This  value  of  ainstr  should  be  consistent 
with  the  standard  deviation  value  used  in  the  Truth  Model.  The  user  also  has 
the  option  of  processing  the  entire  data  batch  or  a  portion  as  well  as  forcing  the 
estimator  to  run  a  specified  number  of  iterations.  Various  other  switches  in  the 
user-input  script  control  the  types  of  output  such  as  plots  and  tables.  Once  the 
program  is  executed,  built-in  checks  are  performed  to  enure  compatible  range  data, 
observational  points,  and  time  data.  The  estimator  will  recursively  execute  the  least- 
squares  algorithm  portion  of  the  code  until  either  all  convergence  criteria  are  satisfied 
or  the  number  of  maximum  iterations  is  reached.  At  that  point,  the  program  will 
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leave  the  iterative  loop  and  propagate  the  estimated  dynamics  based  on  the  last 
valid  estimated  Xre/(£0)-  Depending  on  the  user’s  output  selection,  various  plots 
and  tables  will  be  displayed  or  printed.  Plots  that  are  available  include:  comparison 
plots  of  actual  and  estimated  dynamics,  actual  and  predicted  range  observations,  and 
the  distribution  of  residuals.  Also  filter  performance  data  is  printed  to  the  screen 
such  as  the  standard  deviation  of  the  residuals,  number  of  iterations,  convergence 
status,  and  covariance  values. 

3.5.5  Convergence  Criteria.  In  the  estimation  program,  two  factors  are 
used  to  judge  convergence  to  an  acceptable  estimate  of  the  true  solution.  The  first 
factor  is  the  relative  comparison  between  the  state  correction  5X  and  its  associ¬ 
ated  covariance  term  Psxu ■  When  successive  solutions  of  5X  he  within  the  error 
ellipsoid  indicated  by  the  covariance  matrix,  no  further  iterations  are  meaningful  to 
compute  [15].  This  convergence  criterion  coded  into  the  algorithm  is  expressed  as 

|dXi(t0)|  <C  k2\fP&xii  (3.28) 

where  the  absolute  value  of  the  state  correction,  hXj(f0),  must  be  much  less  that 
the  square  root  of  the  diagonal  terms  of  the  covariance  matrix  Psxu  scaled  down  by 
the  coefficient  fc2 .  Typical  values  explored  for  A:2  compare  dX,(t0)  to  1%  to  50%  of 
J Psxtl  •  Specific  ko  values  will  be  stated  in  the  next  chapter  for  actual  cases.  The 
diagonal  values  of  the  covariance,  Psxu ,  effectively  correspond  to  the  ith  state  which 
is  being  compared  to  the  corresponding  ith  state  in  hXj(A0)- 

The  second  parameter,  used  for  ensuring  convergence,  is  the  statistical  distri¬ 
bution  of  the  residuals  from  the  estimated  trajectory.  Residuals  that  are  sufficiently 
small  enough  will  result  in  a  Gaussian  distribution  that  has  a  standard  deviation 
comparable  to  the  original  standard  deviation  of  the  data  batch.  A  conditional 
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statement  is  coded  in  the  algorithm  to  check  if  the  following  relation  is  true 


^instr  ~  &resid  yo.Zuj 

where  ainstr  is  the  standard  deviation  associated  with  the  instrument  measurements 
defined  in  the  Q  matrix.  The  standard  deviation  of  the  residuals  from  the  estimation 
is  (Jresid ■  The  A'3  scaling  coefficient  is  used  to  specify  a  percentage  of  error  that 
is  acceptable  for  Equation  3.29.  Although,  <5Xj(A0)  may  be  relatively  small  and 
meet  the  first  criterion  above,  the  residuals  may  indicate  a  different  story  since 
the  residuals  are  not  dependent  on  the  covariance  P§x-  In  a  common  case  where 
the  first  criterion  is  met  but  the  residuals  are  not  sufficiently  small,  results  in  a 
<jresid  that  is  near  or  approaching  crmstr .  This  condition  is  easily  identified  by  a 
Gaussian  distribution  of  the  residuals  which  appears  to  have  a  non-zero  mean.  The 
Aq  coefficient  is  chosen  such  that  the  mean  of  the  residual  distribution  is  near  zero 
satisfying  Equation  3.29  above.  Allowing  the  algorithm  to  perform  a  couple  more 
successive  iterations  typically  corrects  this  non-zero  mean  problem  and  results  in  a 
fully  converged  solution. 

Stability  of  the  estimator  is  controlled  by  the  Aq  scaling  coefficient.  This  coef¬ 
ficient  scales  the  state  correction  variables  proportionally  by  a  constant.  A  modified 
form  of  Equation  2.4.3  is  used  in  the  algorithm  as  shown  below 

SX(t0)  =  hPsx  ^Q~lvt  (3.30) 

i 

Cases  exist  which  the  estimator  fails  to  converge  on  a  solution  by  inadvertently 
maximizing  the  S  function,  easily  identified  by  residuals  approaching  infinity.  By 
proportionally  scaling  down  the  state  correction  vector  using  Aq ,  the  algorithm  be¬ 
haves  more  stable  by  making  smaller  corrections  and  typically  moves  toward  the 
desired  solution,  albeit  through  more  iterations  than  would  normally  be  required. 
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IV.  Simulation  &  Results 


Beginning  with  a  hypothetical  scenario  that  provides  a  basis  for  the  a  priori  knowl¬ 
edge  needed  for  this  estimation  process,  this  Chapter  will  explore  three  estimation 
cases  through  simulation.  The  initial  conditions  formed  from  the  scenario  will  be 
the  basis  for  the  estimation  cases.  The  three  cases  will  explore  the  effects  of  the 
following  on  the  success  of  the  nonlinear  least  squares  estimator:  noisy  data,  the  a 
priori  estimates,  and  the  number  of  observational  points  used  (i.e.  number  of  data 
batches  processed).  Each  case  will  result  in  estimates  of  the  rotational  dynamics 
and  moments  of  inertia,  which  will  provide  insight  into  the  dynamical  and  physical 
properties  of  the  hypothetical  resident  space  object. 

4-1  Scenario 

Consistent  with  the  assumptions  made  in  the  derivation  of  the  equations  of 
motion  for  the  system  dynamics,  the  Target  is  assumed  to  be  in  a  true  circular  orbit, 
free  from  external  perturbations,  for  the  duration  of  the  simulated  proximity  opera¬ 
tion  in  which  observational  data  is  being  collected.  The  Target  is  said  to  be  in  LEO 
with  an  altitude  of  800  km  (497 mi).  In  nonlinear  least  squares  estimations,  having 
initial  estimates  or  a  priori  knowledge  of  the  state  parameters  is  required.  Choosing 
the  initial  values  for  these  parameters  is  absolutely  critical  in  the  estimation  outcome 
(i.e.  convergence  vs.  divergence)  as  well  as  the  quality  of  the  estimated  solution  (i.e. 
believable  or  physically  impossible).  Although  it  is  the  goal  to  ultimately  estimate 
the  rotational  dynamics  and  moments  of  inertia  of  the  RSO,  the  estimation  process 
cannot  start  without  some  initial  information  for  these  parameters.  Ideally,  if  the 
target  is  categorized  as  cooperative ,  the  technical  data  of  the  RSO  would  be  consid¬ 
ered  available  and  known  with  relatively  high  accuracy.  However,  for  uncooperative 
RSO  targets,  this  may  pose  a  challenge  requiring  some  creativity  for  developing  the 
initial  estimates  and  associated  physical  information  (i.e.  mass,  moments  of  inertia, 
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dimensions,  structure,  etc).  In  such  cases  several  methods  and  tools  exist  to  aid  in 
the  development  of  this  information.  One  method  is  space  surveillance  using  ground 
radar.  FGAN  Research  Institute  for  High-Frequency  Physics  and  Radar  Techniques 
in  Wachtberg,  Germany,  utilizes  a  34-meter  Tracking  and  Imaging  Radar  (TIRA)  to 
investigate  radar  techniques  for  space  surveillance  applications  [26] .  One  application 
is  providing  attitude  and  configuration  information  for  anomaly  resolution  [26] .  After 
the  Japanese  Advanced  Earth  Observation  Satellite-I  (ADEOS-I)  spacecraft  experi¬ 
enced  a  malfunction  in  1997,  FGAN  was  able  to  use  radar  imagery  of  the  spacecraft 
along  with  a  wire- grid  model  analysis  to  determine  the  loss  of  the  0.5  mm-thick  solar 
panel  as  well  as  the  spacecrafts  angular  velocity  components  [26,  27,  28]. 


Figure  4.1  Radar  Image  of  the  ADEOS  I  Spacecraft  (courtesy  of  FGAN-FHR) 


FGAN  determined  that  the  spacecraft  was  rotating  about  two  axes  with  angu¬ 
lar  velocities  of  approximately  0.1  deg/ s  and  0.4  deg/ s  [26,  28].  Using  these  quantities 
solely  as  a  priori  knowledge  for  this  scenario,  the  angular  velocities  will  be  arbitrarily 
assigned  to  the  following  axes  in  the  body  frame  of  the  Target. 


rad  ,  „  rad 

uji  =  0.0069 -  and  nq  =  0.0017 - 


(4.1) 


s 
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Strictly  for  visualization  purposes,  the  ADEOS  I  spacecraft  will  represent  the 
generic  RSO  Target  body  used  in  this  simulation.  Its  distinctive  geometry  and 
partially  known  characteristics  make  it  well  suited  for  this  scenario.  Therefore,  Fig¬ 
ure  4.2  illustrates  the  Target  along  with  the  assigned  body  axes  in  §7,  and  the  obser¬ 
vational  points  (n  =  1,  2,  3)  added. 


bi 


Figure  4.2  Simulated  Resident  Space  Object 

With  the  orbital  information,  angular  velocity  information,  observation  point 
locations,  and  the  Target  geometry  illustrated,  the  states  at  epoch  are  now  pro¬ 
vided.  The  rotational  states  of  the  Target  at  time  t0  are  shown  in  Table  4.1;  for 
simplicity,  the  body  frame  and  the  orbital  frame  are  aligned  at  epoch.  The  relative 


Table  4.1  Rotational  States  of  the  RSO  at  Epoch 


[rad/s] 

u>2  [rad/s] 

[rad/s] 

9i  [rad] 

02  [rad] 

03  [rad] 

A  Priori 

0.0069 

0 

0.0017 

0 

0 

0 

Truth 

0.01 

0 

0.001 

0 

0 

0 

position  between  the  Target  and  the  Observer  Spacecraft  should  be  relatively  small 
for  a  proximity  operation.  With  the  example  of  the  RELAVIS  LIDAR  sensor  which 
has  a  maximum  measurement  range  of  5  km  [12]  in  mind,  the  relative  motion  states 
are  provided  in  Table  4.2.  Although  this  particular  Target  would  realistically  be 
described  as  semi-rigid,  for  the  purpose  of  this  thesis  and  consistent  with  the  de¬ 
terministic  dynamics  developed  in  Chapter  2,  the  RSO  Target  is  considered  a  rigid 
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Table  4.2  Relative  Position  &  Velocity  States  of  the  RSO  at  Epoch 


Sx  [m\ 

Sy  [m] 

Sz  [m] 

Sx  [m/s] 

Sy  [m/s] 

Sz  [m/s] 

A  Priori 

54 

55 

49 

0.01 

0.01 

0.01 

Truth 

50 

50 

50 

0.01 

0.01 

0.01 

Error 

8% 

10% 

2% 

0% 

0% 

0% 

body  for  the  duration  of  the  observation,  having  constant  moments  of  inertia  shown 
in  Table  4.3.  Three  trials  will  be  accomplished  for  the  9-state  estimation  cases.  Each 


Table  4.3  Moments  of  Inertia  of  RSO 


A  [ kgm 2] 

B  [kgm2] 

C  [kgm2] 

Error 

Trial  1 

50.5 

33 

13.5 

10% 

Trial  2 

54 

36 

12 

20% 

Trial  3 

58.5 

39 

10.5 

30% 

Truth 

45 

30 

15 

- 

trial  will  use  initial  conditions  that  consistently  vary  in  percent  error  from  the  true 
moment  of  inertia  value. 

The  tracked  observational  points  used  on  the  Target  body,  measured  from 
the  assumed  center-of-mass  location,  are  illustrated  in  Figure  4.2  and  provided  in 
Table  4.4.  Since  the  observation  points  are  not  being  estimated,  no  distinction  is 
made  between  their  true  location  and  their  approximated  location. 


Table  4.4  Observational  Point  Locations  on  RSO  Body 


n 

bln  [m] 

b2n  [m] 

hn  M 

i 

0 

0 

20 

2 

2 

6 

-1 

3 

4 

-1 

-1 

With  the  true  initial  conditions  provided  in  the  tables  above,  the  Truth  Model 
program  is  executed  resulting  in  the  true  dynamics  model  of  the  Observer- Target 
system  as  well  as  generated  data  batches  for  simulated  range  measurements.  The 
relative  position  of  the  Observer  with  respect  to  the  Target  is  illustrated  in  Figure  4.3. 
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Figure  4.4  depicts  the  relative  motion  of  the  Observer  about  the  Target  for  one  orbital 
period  (of  the  Target).  Figure  4.5  shows  the  trace  that  the  tracked  observational 


Observer  Spacecraft  Position  Relative  to  Target  [e-frame] 
[1  orbital  period=100.8667min] 


Figure  4.3  Relative  Position  Multi-View 

points  make,  viewed  in  $e,  as  the  Target  body  rotates  with  angular  velocities,  uj\ 
and  u>3,  during  the  span  of  one  orbital  period. 
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Cross-Track  [m] 


Spacecraft  Motion  Relative  to  Target 
(100.8667min  from  Epoch) 


Radial-Track  [m],  e^ 

Figure  4.4  Observer  Spacecraft  &  RSO  Relative  Motion 


Trace  of  Tracked  Observational  Point(s)  on  Rotating  Target  Body 
(100.8667min  from  Epoch) 


Figure  4.5  Observational  Point  Trace  of  RSO  Body 
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4-2  Case  I 


This  case  begins  with  the  use  of  the  RSS  surface  plot  method  to  illustrate  the 
effects  of  the  number  of  points  (i.e.  data  sets)  on  local  minima  and  solution  ambi¬ 
guity.  Following  the  RSS  analysis  and  generated  best-guess  values,  one  estimation 
attempt  is  made  using  data  with  no  noise.  This  case  involves  the  following: 

•  RSS  Surface  Plot  Analysis 

•  Data  Type:  Perfect  (no  noise) 

•  Number  of  Observational  Points:  One  for  Estimation 

•  Observational  Time  Span:  A  Orbital  Period  of  Target 

•  Estimated  States:  12-State  Version 

4-2.1  RSS  Analysis.  A  method  of  refining  an  initial  estimate  is  the  use  of 
a  residual  sum  squared  surface  plot.  In  the  RSS  program,  the  surface  plot  provides 
an  RSS  value  as  a  function  of  two  state  variables.  The  RSS  value  is  essentially  a 
measure  of  how  well  the  estimate  fits  the  actual  data.  The  two  parameter  values 
that  provide  the  minimum  RSS  value  are  the  initial  guess  values  that  should  be  used 
in  the  estimation  algorithm. 

Consider  the  that  initial  estimates  for  the  relative  position  and  velocity  states 
are  representative  of  their  true  values  to  some  satisfactory  level  of  accuracy.  However, 
also  consider  that  the  estimates  of  the  rotational  velocities,  determined  by  radar 
image  analysis,  is  considered  a  low-confidence  estimate.  Using  the  RSS  analysis, 
the  low-confidence  estimate  can  be  used  to  develop  an  estimate  that  is  closer  to 
the  actual  value.  Table  4.1  provides  the  rough  estimates  based  on  the  radar  image 
analysis.  With  the  accuracy  of  the  analysis  known,  the  range  or  bound  for  the  guess 
is  defined.  For  the  purpose  of  this  example,  consider  that  the  image  analysis  values 
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lead  one  to  believe  the  true  values  lie  in  the  following  range 


0.001  <  fjj\  (to)  <  0.02  and  —  0.0005  <  uq (to)  <  0.05  (4.2) 

and  define  uq  as  parameter  1  and  u 3  as  parameter  2.  The  user’s  inputs  include  the 
a  priori  estimates  for  relative  motion  from  Table  4.2,  the  observational  points  from 
Table  4.4,  and  the  settings  shown  in  Table  4.5.  Since  running  this  program  with  these 


Table  4.5  RSS  Guess  Generator  Boundary  Conditions 


Parameter 

Lower  Bound 

Upper  Bound 

No.  of  Increments 

1 

0.001 

0.020 

100 

2 

-0.0005 

0.050 

100 

settings  for  all  data  sets  takes  approximately  eight  to  ten  hours  of  CPU  time,  and 
with  memory  limitations  of  the  computer  used,  only  one  tenth  of  the  total  data  batch 
(606  data  points)  will  be  used  as  the  sample  size  for  the  RSS  analysis.  This  provides 
an  adequate  sample  for  determining  the  initial  conditions  of  the  two  parameters  for 
this  case.  The  surface  plots  and  their  associated  contour  plots  resulting  from  the 
execution  of  this  program  are  shown  in  the  Figure  4.6. 

Figures  4.6a  and  4.6b  show  the  computed  results  of  the  first  RSS  solution 
(Run  1)  using  observation  point  number  one  (i.e.  one  data  batch).  Recall  that  the 
estimator  will  strive  to  minimize  the  sum  of  the  residuals  squared  (i.e.  the  RSS 
value)  by  finding  a  local  minimum.  Within  the  guess  bounds  (Equ.  4.2),  two  distinct 
local  minimums  exist  as  clearly  seen  in  the  surface  plot.  One  can  see  from  the 
contour  plot  (Fig.  4.6b)  that  the  a  priori  estimates  for  uq  and  u>3  (Table  4.1)  lies 
roughly  in  between  the  minimums  (the  dark  regions).  This  may  pose  a  problem 
for  the  estimator,  as  it  may  choose  the  wrong  minimum,  where  the  desired  solution 
does  not  he.  Keep  in  mind  the  location  of  the  true  solution  as  annotated  on  the 
contour  figure  and  given  in  Table  4.1.  The  global  minimum  is  determined  by  the 
RSS  program  to  be  the  point  where  parameter  1  (uq)  is  0.0050  and  parameter  2  (cu3) 
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is  0.0102.  This  can  be  seen  as  the  region  which  dips  the  lowest  in  surface  plot  and 
the  other  major  dark  region  in  the  contour  plot.  However,  this  is  deceivingly  the 
wrong  minimum  and  may  very  well  lead  the  estimator  in  the  wrong  direction  to  an 
incorrect  solution. 

By  re-executing  the  RSS  program  (Run  2)  using  observation  points  one  and 
two  (i.e.  two  independent  data  batches),  the  program  provides  Figures  4.6c  and  4.6d. 
Notice  that  the  global  minimum  from  the  first  run  (Figure  4.6a  and  4.6b)  is  being 
reduced  in  depth  and  that  the  second  local  minimum  (where  the  solution  actually 
lies)  is  now  the  global  minimum  in  this  run. 

Run  3  of  the  RSS  program,  using  a  total  of  three  observational  points,  the 
program  results  in  the  plots  seen  in  Figures  4.6e  and  4.6f.  Notice  that  in  Figure  4.6e 
only  one  major  dip  (minimum)  is  seen.  From  a  contour  perspective,  this  is  clearly 
seen  in  Figure  4.6f,  as  only  one  major  and  defined  dark  region  exists  in  the  location 
of  the  true  solution.  This  region  is  where  the  program- determined  guess  lies,  which 
is  extremely  close  to  the  true  location  of  the  solution. 

A  summary  of  the  results  from  the  three  RSS  program  runs  are  provided  in 
Table  4.6.  The  results  are  shown  up  to  eight  decimal  digits  to  illustrate  the  imper¬ 
fection.  The  number  of  observation  points  used  correspond  to  the  points  listed  in 
Table  4.4  starting  with  n  =  1.  The  percent  error  between  the  generated  guess  for 
ui  and  u 3  and  the  true  values  (Table  4.1)  are  also  given.  By  introducing  multiple 
independent  batches  of  data  from  multiple  range  observations,  the  RSS  solution  is 
drastically  improved. 


Table  4.6  RSS  Parameter  Results 


Run  No. 

Obs  Points  Used,  n 

Par.  1  Estimate 

Error 

Par.  2  Estimate 

Error 

1 

1 

0.00503030 

50% 

0.01021212 

921% 

2 

1,2 

0.01040404 

4% 

0.00256061 

156% 

3 

1,2,3 

0.01002020 

0.2% 

0.00103030 

3% 
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RSS  Solution  Surface  Plot 


Contour  Plot  of  RSS  Solution 


(a) 


RSS  Solution  Surface  Plot 


-0.02  0 


(c) 


RSS  Solution  Surface  Plot 


0.002  0.004  0.006  0.008  0.01  0.012  0.014  0.016  0.018 
Parameter  1 


(b) 


Contour  Plot  of  RSS  Solution 


0.002  0.004  0.006  0.008  0.01  0.012  0.014  0.016  0.018 
Parameter  1 


(d) 


0.002  0.004  0.006  0.008  0.01  0.012  0.014  0.016  0.018 
Parameter  1 


(e)  (f) 

Figure  4.6  RSS  Surface  &  Contour  Plots 
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Theoretically,  the  program  should  give  a  solution  for  two  parameters  (state 
variables)  that  result  in  an  RSS  value  of  zero,  a  perfect  data  fit.  However,  recall 
from  Chapter  3  that  this  program  is  limited  to  providing  a  solution  to  two  state 
variables.  In  this  case  because  there  are  12  total  states,  10  of  the  variables  are 
considered  known  constants  by  the  RSS  program.  The  imperfect  results  seen  in 
Table  4.6  are  a  consequence  of  imperfect  initial  values  for  those  other  10  states  used 
in  generating  this  result.  The  results  from  this  RSS  analysis  will  serve  as  the  initial 
state  estimates  for  the  estimation  runs  to  follow. 

The  estimation  attempt  in  this  case  will  make  use  of  the  results  from  Run  3. 
Table  4.7  provides  the  values  used  in  Case  I. 


Table  4.7  Case  I:  Refined  Rotational  States  of  the  RSO  at  Epoch 


uj\  [rad/s] 

(jl>2  [rad/s] 

u 3  [rad/s] 

^3 

Q 

.  K 

^8 

02  [rad] 

O3  [rad] 

Case  1A 

0.0100 

0 

0.0010 

0 

0 

0 

Truth 

0.01 

0 

0.001 

0 

0 

0 

Notice  Run  3  of  the  RSS  program  (Table  4.6)  produced  values  for  uq  and 
u3  that  are  essentially  equal  to  the  true  state  values  up  to  three  significant  digits. 
In  order  to  fully  make  the  estimation  cases  worthwhile  and  to  legitimately  test  the 
estimator’s  performance,  values  near  the  true  solution  but  not  equal  to,  are  preferred 
in  this  study.  Therefore,  Table  4.8  lists  the  refined  a  priori  initial  conditions  (from 
Run  2)  used  for  the  Case  II  and  Case  III  estimations. 


Table  4.8  Refined  Rotational  States  at  Epoch  for  Cases  II  and  III 


uj\  [rad/s] 

UO2  [rad/s] 

CJ3  [rad/s] 

61  [rad] 

O2  [rad] 

03  [rad] 

A  Priori 

0.0104 

0 

0.0026 

0 

0 

0 

Truth 

0.01 

0 

0.001 

0 

0 

0 

4-2.2  Data  Batch.  The  data  set  processed  in  this  case  uses  a  606-observation 
sample  of  the  entire  data  batch,  or  one-tenth  of  the  Target’s  orbital  period  worth. 
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This  is  done  only  to  be  consistent  with  the  RSS  Analysis,  which  also  used  the  first 
606  measurements.  As  mentioned,  this  case  will  make  use  of  perfect  simulated  data 
containing  no  Gaussian  noise.  The  point  number  used  corresponds  to  the  observa¬ 
tional  points  listed  in  Table  4.4.  Table  4.9  provides  a  summary  of  the  data  set  used 
in  this  case. 


Table  4.9  Case  I:  Data  Batch  Summary 


Point  No.  n 

Batch  Size 

Start  Time  [s] 

End  Time  [s] 

Time  Step  [ s ] 

i 

606 

0 

605 

i 

4-2.3  Estimation  Case  1A.  This  portion  of  Case  I  uses  the  best-guess  so¬ 
lution  for  u>i  and  cu3  provided  by  Run  3  of  the  RSS  surface  plot  analysis.  The  initial 
rotational  dynamics  states  are  provided  in  Table  4.7.  Using  the  initial  conditions  for 
relative  position  and  velocity  listed  in  Table  4.2  and  observation  point  number  one 
listed  in  Table  4.4,  the  estimator  program  is  executed  with  the  following  setting  pro¬ 
vided  in  Table  4.10  below.  The  setting  imax  represents  the  user-specified  maximum 


Table  4.10  Case  I:  Estimator  Settings  for  12-State  Version 


^ max 

&instr 

ki 

fa 

fa 

States 

NumPts 

NoiseMode 

200 

5  x  10“5 

0.25 

0.01 

0.20 

12 

i 

off 

number  of  iterations  allowed  before  the  program  automatically  terminates,  States 
defines  the  mode  (9  or  12)  for  the  estimated  state  version  (Equ.  2.65  or  Equ.  2.70), 
NumPts  sets  the  number  of  points  to  be  used  in  the  estimation  process  (i.e.  range 
data  batches),  and  the  NoiseMode  string  simply  has  the  program  use  the  version 
of  the  range  data  containing  no  noise.  Since  the  data  contains  no  noise,  it  has  a 
standard  deviation  of  <Jinstr  =  0.  However,  to  avoid  the  Q  matrix  (Equ.  2.97)  from 
becoming  singular,  <Jjnstr  =  0.00005  is  used  instead. 

The  estimator  program  converged  on  a  solution  in  38  iterations.  Based  on  the 


predicted  range  measurements,  no  significant  difference  exists  between  the  estima- 
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tor’s  solution  and  the  actual  solution.  Table  4.11  compares  the  true  values  of  the 


Actual  Range  Measurements  vs.  Predicted  Range  Measurements 
Target  Observation  Point  1 


Figure  4.7  Case  I:  Predicted  Range  Profile  (12-State) 


states  at  epoch  and  the  estimated  results  for  this  case.  Propagating  the  estimator’s 


Table  4.11  Case  I:  Estimated  12-State  Solution  at  Epoch 


State  Variable 

Estimate,  X(£q) 

Truth,  X(£o) 

Error 

5x 

[m\ 

52.16 

50 

4.3% 

8y 

[m\ 

51.84 

50 

3.7% 

5z 

M 

46.42 

50 

7.2% 

5x 

[m/s] 

0.0163 

0.01 

63% 

Sy 

[m/s] 

-0.0025 

0.01 

125% 

8z 

[m/s] 

0.0152 

0.01 

52% 

LOl 

[rad/ s] 

0.0101 

0.01 

0.0001 

<^2 

[rad/ s] 

0.0001 

0 

0.0001 

^3 

[rad/ s] 

0.0010 

0.001 

0.0001 

[rad] 

-0.1167 

0 

0.1167 

02 

[rad] 

-0.0224 

0 

0.0224 

@3 

[rad] 

-0.0830 

0 

0.0830 

solution  for  the  states  at  epoch,  results  in  the  dynamics  shown  for  relative  position 

in  Figure  4.8.  Figure  4.9  illustrates  the  Euler  angles.  Additional  data  is  presented 
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Time  [s] 


Figure  4.8  Case  I:  Estimated  Relative  Position  States  (12-State) 


in  Appendix  C  for  this  case.  It  should  be  mentioned  that  in  this  estimation  run, 


Euler  Angle  Orientation  of  Target 


Time  [s] 

Figure  4.9  Case  I:  Estimated  Euler  Orientation  Angle  States  (12-State) 


had  the  aresid  setting  been  changed  to  a  value  closer  to  machine  zero,  this  run  would 
have  failed  the  standard  deviation  convergence  criterion.  This  solution  is  less  ac¬ 
curate  than  ideally  desired  for  a  case  without  noise  and  relatively  accurate  initial 
values.  The  major  error  ranging  from  52%  to  125%  is  seen  in  the  relative  velocity 
states  (Table  4.11).  These  errors  and  <Jresid  are  due  to  the  initial  conditions  of  the 
relative  position  variables.  This  estimation  run  is  considered  successful  and  satisfied 
its  purpose  to  illustrated  a  run  using  the  RSS  analysis  results  for  data  processing. 
Further  cases  in  this  research  will  make  use  of  the  entire  data  batch  containing  noise. 
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4-3  Case  II 


This  will  be  the  first  case  in  which  the  estimation  process  is  applied  to  simulated 
data  containing  Gaussian  noise.  This  case  will  involve  the  following: 

•  Data  Type:  Noisy  (with  Gaussian  Error) 

•  Number  of  Observational  Points:  One 

•  Observational  Time  Span:  1  Orbital  Period 

•  Estimated  States:  12-State  &  9-State  Versions 

The  objective  of  this  case  is  to  generate  an  estimated  solution  using  a  single  tracked 
observation  point  (i.e.  one  range  data  set)  for  comparison  with  Case  III  involving 
multiple  observation  points.  Additional  data  is  provided  in  Appendix  C. 

4-3.1  Data  Batch.  The  generated  data  batch  that  is  processed  in  both 
portions  of  this  case  has  6053  simulated  range  measurements  with  a  standard  de¬ 
viation  listed  below.  A  negligible  difference  between  the  standard  deviation  of  the 
instrument  and  the  actual  standard  deviation  of  the  data  exists.  This  is  simply 
due  to  the  user’s  input  of  <Jinstr  in  the  Truth  Model  script  and  the  pseudo-random 
number  algorithm  generating  data  nearly  equal  with  some  difference  because  of  the 
randomness  involved.  The  value  of  crinstr  will  be  used  in  estimation  process.  The 
point  number  used  corresponds  to  the  observational  points  listed  in  Table  4.4. 


Table  4.12  Case  II:  Data  Batch  Summary 


Point  No.  n 

®instr 

& act 

Batch  Size 

Start  Time  [s] 

End  Time  [ s ] 

Time  Step  [s] 

1 

0.250 

0.248 

6053 

0 

6052 

1 

4-3.2  12- State  Estimation.  This  portion  of  the  case  will  involve  estimating 

both  the  relative  motion  dynamics  and  the  rotational  states  of  the  RSO.  Recall  the 
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12-state  vector  defined  as, 


X 


Sx  Sy  Sz  Sx  Sy  Sz  u>i  u>2  <^3  #1  #2 


03 


T 


During  several  initial  runs  made  during  the  development  of  the  estimator  code,  it 
was  determined  that  the  program  is  generally  stable  and  performs  the  best  when 
0.01  <  k\  <  0.50.  In  this  case,  the  estimator  setting  are  provided  below 


Table  4.13  Case  II:  Estimator  Settings  for  12-State  Version 


T'max 

&instr 

ki 

States 

NumPts 

NoiseMode 

200 

0.25 

0.50 

0.05 

0.012 

12 

1 

on 

With  the  settings  above,  the  program  ran  for  approximately  4  hrs  where  it 
reached  the  maximum  allowed  number  of  iterations  without  convergence  to  an  ac¬ 
ceptable  solution.  The  covariance  criterion  was  met;  however,  the  standard  devia¬ 
tion  criterion  was  not.  The  standard  deviation  of  the  residuals  from  the  estimation 
process  was  approximately  crresid  ~  14,  as  seen  in  Figure  4.15  while  the  standard 
deviation  of  noise  error  in  the  data  batch  is  (Tinfstr  ~  0.25.  Therefore,  this  refer¬ 
ence  solution  remains  only  a  reference  solution  and  is  not  considered  the  estimated 
trajectory.  The  last  reference  solution,  Xre/(f0),  prior  to  program  termination  is 
provided  in  the  Table  4.14  along  with  the  associated  percent  relative  error  from  the 
true  states.  To  illustrate  what  the  non-convergent  solution  resembles,  the  resulting 
reference  solution  at  imax  is  propagated  forward  in  time  to  the  end  of  the  observa¬ 
tional  period.  The  relative  position  and  velocity  histories  are  shown  in  Figures  4.10 
and  4.11.  Their  respective  true  error  as  a  function  of  time  is  shown  in  Figures  C.8 
and  C.9  in  Appendix  C.  The  six  rotational  dynamical  states  of  the  Target  are  il¬ 
lustrated  in  Figures  4.12  and  4.13.  Their  associated  true  error  plots  are  shown  in 
Figures  C.10  and  C.ll.  The  overall  solution  of  the  reference  trajectory  resulted 
in  the  predicted  range  observations  shown  in  Figure  4.14  compared  with  the  actual 

noisy  range  data  set.  Notice  that  the  solution  results  in  a  predicted  range  profile 
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Figure  4.10  Case  II:  Estimated  Relative  Position  States  (12-State) 


Relative  Spacecraft  Velocity  WRT  Target 


Figure  4.11  Case  II:  Estimated  Relative  Velocity  States  (12-State) 
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Figure  4.12  Case  II:  Estimated  Angular  Velocity  States  (12-State) 
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Table  4.14  Case  II:  Estimated  12-State  Results  at  Epoch 


State  Variable 

Estimate,  Xrey(£o) 

Truth,  X(£o) 

Error 

5x 

[to] 

32.06 

50 

35.9% 

Sy 

[to] 

-52.43 

50 

204.9% 

5z 

[m\ 

76.27 

50 

52.5% 

5x 

[m/s] 

-0.0022 

0.01 

122% 

8y 

[m/s] 

0.0409 

0.01 

309% 

Sz 

[m/s] 

0.0159 

0.01 

59% 

CJl 

[rad/ s] 

0.0248 

0.01 

0.0148 

UJ2 

[rad/ s] 

-0.0434 

0 

0.0434 

ids 

[rad/ s] 

-0.0402 

0.001 

0.0392 

<h 

[rad] 

4.7289 

0 

4.7289 

@2 

[rad] 

0.9465 

0 

0.9465 

@3 

[rad] 

8.5974 

0 

8.5974 

Figure  4.13  Case  II:  Estimated  Euler  Orientation  Angle  States  (12-State) 


that  does  generally  resemble  the  true  range  observations.  The  true  range  profile  is 
more  nonlinear  than  it  appears  in  Figure  4.14  simply  because  of  the  plot  scale.  The 
estimated  solution  fits  the  middle  portion  without  much  trouble;  however,  the  most 
nonlinear  regions,  located  approximately  from  0  —  2000  s  and  4000  —  6050  s,  posed 
the  greatest  difficulty  for  the  estimator  in  this  case.  The  distribution  of  the  residu¬ 
als  between  the  predicted  range  observations  and  the  actual  observations  provides  a 
means  to  judge  the  quality  of  the  solution.  Since  the  noise  or  error  of  the  simulated 
measurements  is  Gaussian  in  nature  (by  creation),  the  residuals  of  the  predicted 
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Actual  Range  Measurements  vs.  Predicted  Range  Measurements 
Target  Observation  Point  1 


Figure  4.14  Case  II:  Predicted  Range  Profile  (12-State) 

observations  should  also  resemble  a  Gaussian  distribution  with  a  similar  standard 
deviation.  However,  in  this  case  crresui  ^  crinstr >  which  is  why  the  convergence  criteria 
failed.  With  the  k 3  setting  defined  in  Table  4.13,  the  standard  deviation  criterion 
required  (Jresid  to  fall  in  the  range  of  0.247  <  <Jresid  <  0.253. 

Figure  4.15a  depicts  a  residual  distribution  that  loosely  resembles  a  distorted 
Gaussian  curve  centered  on  zero,  but  with  a  standard  deviation  of  approximately  14 
meters.  In  Figure  4.15b  the  distribution  of  range  residuals  as  a  function  of  observa¬ 
tion  time  is  shown  to  have  a  desirable  mean  near  zero,  as  mentioned;  however,  the 
spread  of  residuals  is  not  characteristic  of  a  true  normal  distribution.  From  iteration 
28  to  200,  the  standard  deviation  of  the  residuals  remained  at  aresid  ~  14.042.  It 
would  appear  that  the  estimator  reached  a  local  minimum  indicated  by  this  appar¬ 
ent  limit  of  the  standard  deviation.  This  will  be  discussed  further  in  the  last  section 
of  this  chapter.  The  covariance  criterion  was  satisfied  at  iteration  55  till  program 
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(a)  (b) 

Figure  4.15  Case  II:  Residual  Distribution  for  Point  No.  1  (12-State) 

termination.  This  estimation  attempt,  that  did  not  converge  in  200  or  less  iterations 
with  the  given  filter  settings,  is  considered  a  failed  attempt. 


4-3.3  9- State  Estimation.  This  portion  of  Case  II  will  attempt  to  estimate 

the  rotational  states  as  well  as  the  moments  of  inertia  of  the  Target.  To  further 
explore  the  ability  to  estimate  the  moments  of  inertia,  three  trials  will  be  undertaken. 
Recall  the  9-state  vector  is  given  as 


X  = 


1  T 


ail  uj2  eo3  9\  62  9%  A  B  C 


The  estimator  settings  for  this  run  are  shown  in  Table  4.15.  This  case  will  make 
Table  4.15  Case  II:  Estimator  Settings  for  9-State  Version 


^ max 

&instr 

h 

k2 

k3 

States 

NumPts 

NoiseMode 

200 

0.25 

0.50 

0.10 

0.012 

9 

1 

on 

use  of  the  same  range  data  batch  from  the  same  single  observation  point  used  in 

the  12-state  run.  The  a  priori  initial  estimates  for  the  rotational  states  are  listed 
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in  Table  4.8.  Unlike  the  12-State  portion,  this  9-State  version  of  this  process  will 
be  completed  for  three  trials  using  the  initial  conditions  listed  in  Table  4.3  for  the 
moments  of  inertia.  Note  that  each  trial  increases  the  relative  error  between  the 
initial  guess  and  the  true  value  of  the  MOI.  After  running  the  estimator  program  for 
each  trial,  keeping  all  initial  states  besides  the  moments  of  inertia  the  same,  results 
in  the  following  shown  in  Table  4.16  In  this  case,  all  trials  satisfied  the  convergence 


Table  4.16  Case  II:  Estimated  9-State  Results  at  Epoch 


Estimates 

Variable 

Trial  1 

Trial  2 

Trial  3 

Truth 

[10%  MOI  Error] 

[20%  MOI  Error] 

[30%  MOI  Error] 

LJl 

[rad/ s] 

0.0100 

0.0100 

0.0100 

0.01 

U2 

[rad/ s] 

0.0000 

0.0000 

0.0000 

0 

CJ3 

[rad/ s] 

0.0000 

0.0000 

0.0000 

0.001 

di 

[rad] 

0.0002 

0.0004 

0.0007 

0 

e2 

[rad] 

0.0013 

0.0015 

0.0000 

0 

[rad] 

0.0020 

0.0015 

0.0014 

0 

A 

[m2kg] 

50.76 

53.91 

58.54 

45 

B 

to 

to 

32.73 

36.10 

38.96 

30 

C 

[■ m2kg ] 

13.19 

12.12 

10.43 

15 

®resid  [m] 

0.249 

0.249 

0.249 

0.247 

iterations 

9 

8 

8 

- 

processing  time  [s] 

98 

100 

120 

- 

Angular  Velocity  of  Target 


Figure  4.16  Case  II:  Estimated  Angular  Velocity  States  (9-State) 
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Angular  Velocity  Error 
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Figure  4.17  Case  II:  True  Error  of  Estimated  Angular  Velocities  (9-State) 


Euler  Angle  Orientation  of  Target 


Figure  4.18  Case  II:  Estimated  Euler  Orientation  Angle  States  (9-State) 
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Figure  4.19  Case  II:  True  Error  of  Estimated  Euler  Orientation  Angles  (9-State) 
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Table  4.17  Case  II:  Percent  Error  of  9-State  Estimates 


Error  of  Estimates 

Variable 

Trial  1 

Trial  2 

Trial  3 

[10%  MOI  Error] 

[20%  MOI  Error] 

[30%  MOI  Error] 

[rad/s] 

0 

0 

0 

UJ2  [rad/s] 

0 

0 

0 

uj 3  [rad/s] 

0.001 

0.001 

0.001 

,  K 

i—l 

0.0002 

0.0004 

0.0007 

62  [rad] 

0.0013 

0.0015 

0 

O3  [rad] 

0.0020 

0.0015 

0.0014 

A 

12.8% 

19.8% 

30.1% 

B 

9.1% 

20.3% 

29.9% 

C 

12.1% 

19.2% 

30.5% 

&resid 

0.8% 

0.8% 

0.8% 

criteria.  Propagating  the  estimated  solution  for  Trial  2  forward  in  time,  yields  the 
following  results.  Using  the  Trial  2  solution  as  an  example,  the  predicted  range 
measurements  are  shown  in  Figure  4.20.  To  illustrate  that  the  predicted  range  is  truly 
a  perfect  fit  to  the  actual  data  set  and  well  within  the  noise,  Figure  C.12  provides  a 
zoomed  in  caption  of  the  range  profile  in  the  region  of  5740  —  6050  s.  The  distribution 
of  the  residuals  is  provided  in  Figure  4.21.  The  histogram  (Fig.  4.21a)  shows  the 
estimation  residuals  having  a  Gaussian  distribution  comparable  to  the  distribution  of 
the  induced  measurement  error  in  the  data  set  (e.g.  Fig.  3.3).  Figure  4.21b  illustrates 
the  distribution  of  the  residuals  as  a  function  of  time  with  crresid  annotated.  This  is 
characteristic  of  the  added  Gaussian  noise  (e.g.  Fig.  3.4).  This  9-state  estimation 
run  is  considered  a  success  using  only  one  data  set  in  the  estimation  process.  The 
only  area  of  concern  is  that  the  solution  indicates  one  axis  of  rotation.  This  will  be 
discussed  in  a  later  section. 
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Actual  Range  Measurements  vs.  Predicted  Range  Measurements 
Target  Observation  Point  1 


Figure  4.20  Case  II:  Predicted  Range  Profile  (9-State) 


(a) 


(b) 


Figure  4.21  Case  II:  Residual  Distribution  for  Point  No.  1  (9-State) 
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4-4  Case  III 

The  12  and  9  state  portions  of  this  case  will  examine  the  results  of  the  esti¬ 
mation  process  using  two  data  sets  from  the  observation  of  two  points  on  the  Target 
body  instead  of  just  one.  Specifically,  this  case  will  involve  the  following: 

•  Data  Type:  Noisy  (with  Gaussian  Error) 

•  Number  of  Observational  Points:  Two 

•  Observational  Time  Span:  1  Orbital  Period 

•  Estimated  States:  12-State  &  9-State  Versions 

The  objective  of  this  case  is  to  determine  if  processing  two  data  batches  (from  two 
observation  points)  affects  the  estimated  solution  versus  the  solution  formed  from 
processing  a  single  data  batch.  This  will  be  applied  to  both  versions  of  the  state 
vector.  Additional  data  is  provided  in  Appendix  C. 

4-4-1  Data  Batch.  Two  data  batches  are  processed  in  this  case.  The  first 
data  batch  associated  with  observation  point  1  is  the  same  data  set  used  in  Case  II. 
The  second  data  set  used  is  associated  with  observation  point  2,  as  listed  in  Table  4.4. 
Table  4.18  provides  the  specifics  for  the  data  sets  used  in  this  case. 


Table  4.18  Case  III:  Data  Batch  Summary 


Point  No.  n 

&instr 

& act 

Batch  Size 

Start  Time  [s] 

End  Time  [s] 

Time  Step  [s] 

i 

0.250 

0.248 

6053 

0 

6052 

i 

2 

0.250 

0.247 

6053 

0 

6052 

i 

4-4-2  12- State  Estimation.  This  portion  of  the  case  will  involve  estimating 

both  the  relative  motion  dynamics  and  the  rotational  states  of  the  RSO  using  two 
sets  of  data.  The  estimator  setting  are  provided  Table  4.19. 
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Table  4.19  Case  III:  Estimator  Settings  for  12-State  Version 


^ max 

&instr 

h 

k2 

k3 

States 

NumPts 

NoiseMode 

200 

0.25 

0.50 

0.05 

0.012 

12 

i 

on 

The  program  ran  for  approximately  100  s  where  it  met  both  the  covariance  and 
standard  deviation  criteria  and  converged  to  a  solution  in  13  iterations.  Table  4.20 
lists  the  results  of  the  estimated  solution  for  the  12  states  at  epoch. 


Table  4.20  Case  III:  Estimated  12-State  Results  at  Epoch 


State  Variable 

Estimate,  X(to) 

Truth,  X(£o) 

Error 

Sx 

[m] 

50.01 

50 

0.02% 

Sy 

[m] 

50.04 

50 

0.08% 

6z 

[m\ 

49.95 

50 

0.1% 

5x 

[m/s] 

0.0100 

0.01 

0 

5y 

[m/s] 

0.0100 

0.01 

0 

Sz 

[m/s] 

0.0100 

0.01 

0 

CJ1 

rad/ s] 

0.0100 

0.01 

0 

CJ2 

rad/ s] 

0.0000 

0 

0 

^>3 

[rad/ s] 

0.0010 

0.001 

0 

0i 

[rad] 

-0.0006 

0 

0.0006 

e2 

[rad] 

0.0004 

0 

0.0004 

o3 

[rad] 

0.0018 

0 

0.0008 

&resid 

n=l  H 

0.249 

0.248 

0.4% 

&resid 

n=2  H 

0.251 

0.247 

1.6% 

The  estimated  states  at  epoch  propagated  forward  in  time  to  the  end  of  the 
observational  period  results  in  Figures  4.22  and  4.23,  along  with  their  true  error 
plots  for  relative  motion  shown  in  Appendix  C.  The  rotational  dynamics  from  the 
estimated  solution  are  shown  in  Figures  4.24  and  4.25.  Applying  the  observation 
function  with  the  estimated  dynamics  results  in  the  predicted  range  observations 
for  each  observational  point  shown  in  Figure  4.26.  The  residuals  from  the  predicted 
range  observation  and  the  actual  range  measurements  have  a  distribution  illustrated 
in  Figures  4.27  and  4.28  for  each  data  batch. 
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Figure  4.22  Case  III:  Estimated  Relative  Position  States  (12-State) 


Relative  Spacecraft  Velocity  WRT  Target 


Figure  4.23  Case  III:  Estimated  Relative  Velocity  States  (12-State) 


Angular  Velocity  of  Target 


Figure  4.24  Case  III:  Estimated  Angular  Velocity  States  (12-State) 
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Euler  Angle  Orientation  of  Target 
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Figure  4.25 


Case  III:  Estimated  Euler  Orientation  Angle  States  (12-State) 


Actual  Range  Measurements  vs.  Predicted  Range  Measurements 
Target  Observation  Point  1 


Actual  Range  Measurements  vs.  Predicted  Range  Measurements 
Target  Observation  Point  2 


(a) 

Figure  4.26 


(b) 


Case  III:  Predicted  Range  Profiles  (12-State) 
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(a)  (b) 

Figure  4.28  Case  III:  Residual  Distribution  for  Point  No.  2  (12-State) 
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4-4-3  9-State  Estimation.  This  portion  of  Case  III  will  attempt  to  form 
an  estimated  solution  of  the  rotational  states  as  well  as  the  moments  of  inertia  of 
the  Target  using  two  data  batches  from  two  observation  points.  This  section  will 
also  perform  three  trials  of  the  estimation  process  for  comparison  with  Case  II.  The 
estimator  settings  for  this  run  are  shown  in  Table  4.21.  This  case  will  use  the  same 


Table  4.21  Case  III:  Estimator  Settings  for  9-State  Version 


^ max 

&instr 

h 

fa 

fa 

States 

NumPts 

NoiseMode 

200 

0.25 

0.50 

0.10 

0.012 

9 

i 

on 

data  batch  from  the  previous  12-state  run.  As  mentioned,  this  case  will  process 
two  data  sets,  one  of  which  is  the  same  data  batch  used  in  Case  II.  The  a  priori 
initial  states  for  the  angular  velocities  and  the  Euler  angles  are  listed  in  Table  4.8. 
The  the  three  trials  are  once  again  accomplished  as  listed  in  Table  4.3.  Running 
the  nonlinear  least  squares  estimator  program  for  each  trial,  resulted  in  the  data 
presented  in  Table  C.l 


Table  4.22  Case  III:  Estimated  9-State  Results  at  Epoch 


Estimates 

Variable 

Trial  1 

Trial  2 

Trial  3 

Truth 

[10%  MOI  Error] 

[20%  MOI  Error] 

[30%  MOI  Error] 

CJl 

rad/ s] 

0.0100 

0.0100 

0.0100 

0.01 

UJ2 

rad/ s] 

0.0000 

0.0000 

0.0000 

0 

Ll>3 

rad/ s] 

0.0010 

0.0010 

0.0000 

.001 

9i 

[rad] 

0.0001 

0.0012 

0.0004 

0 

e2 

[rad] 

-0.0003 

0.0002 

-0.0079 

0 

O3 

[rad] 

0.0009 

-0.0023 

-0.0008 

0 

A 

[m2kg] 

46.96 

52.92 

58.21 

45 

B 

[m2kg] 

38.66 

35.60 

39.90 

30 

C 

[m2kg] 

8.28 

17.30 

9.25 

15 

&resid- 

n=l  H 

0.249 

0.252 

0.256 

0.248 

&resid 

n=2  H 

0.249 

0.248 

0.326 

0.247 

iterations 

18 

21 

200 

- 

processing  time  [s] 

213 

268 

2871 

- 
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Table  4.23  Case  III:  Percent  Error  of  9-State  Estimation  Results 


Percent  Error  of  Estimates 

Variable 

Trial  1 

Trial  2 

Trial  3 

[10%  MOI  Error] 

[20%  MOI  Error] 

[30%  MOI  Error] 

CJl 

[rad/ s] 

0 

0 

0 

U>2 

[rad/ s] 

0 

0 

0 

u;3 

[rad/ s] 

0 

0 

0.0010 

0i 

[rad] 

0.0001 

0.0012 

0.0004 

@2 

[rad] 

0.0003 

0.0002 

0.0079 

e3 

[rad] 

0.0009 

0.0023 

0.0008 

A 

[m2kg] 

4.4% 

17.6% 

29.4% 

B 

[■ m2kg ] 

28.9% 

18.7% 

33% 

C 

[m2kg] 

44.8% 

15.3% 

38.3% 

&residn=i 

0.4% 

1.6% 

3.2% 

(Jresidn=2 

0.8% 

0.4% 

32% 

Using  the  Trial  2  estimated  solution,  the  rotational  dynamics  is  shown  in  Fig¬ 
ures  4.29  and  4.31,  along  with  their  associated  true  error  as  a  function  of  time. 
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Figure  4.29  Case  III:  Estimated  Angular  Velocity  States  (9-State) 


The  estimated  dynamics  results  in  a  predicted  range  measurement  for  each 
observational  point  shown  in  Figure  4.33.  The  associated  residuals  are  illustrated  in 
the  Gaussian  curves  shown  in  Figure  4.34  and  Figure  4.35.  The  estimation  residuals 
and  the  measurement  error  satisfied  the  condition  aresici  ~  ainstr. 
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Angular  Velocity  Error 
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Figure  4.30  Case  III:  True  Error  of  Estimated  Angular  Velocities  (9-State) 
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Figure  4.31  Case  III:  Estimated  Euler  Orientation  Angle  States  (9-State) 
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Figure  4.32  Case  III:  True  Error  of  Estimated  Euler  Orientation  Angles  (9-State) 
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Actual  Range  Measurements  vs.  Predicted  Range  Measurements 
Target  Observation  Point  1 


Actual  Range  Measurements  vs.  Predicted  Range  Measurements 
Target  Observation  Point  2 


(a)  (b) 

Figure  4.33  Case  III:  Predicted  Range  Profiles  (9-State) 
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Figure  4.34  Case  III:  Residual  Distribution  for  Point  No.  1  (9-State) 
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(a)  (b) 

Figure  4.35  Case  III:  Residual  Distribution  for  Point  No.  2  (9-State) 

4-5  Discussion 

With  the  results  of  the  simulation  and  estimation  cases  presented,  this  section 
will  provide  a  comprehensive  discussion  of  the  results  through  the  comparison  of 
solution  error  and  estimator  performance  parameters. 

4-5.1  Case  I.  The  RSS  surface  and  contour  plots  clearly  show  how  intro¬ 
ducing  multiple  data  batches  reduces  the  number  of  overall  solutions  that  satisfy 
the  observation  relation  (Equ.  3.11).  Recall  that  Equ.  3.11  has  three  squared  terms, 
which  allows  for  multiple  solutions  to  exist.  Table  4.6  shows  that  for  one  data  set  the 
error  between  the  estimated  and  true  solutions  for  parameters  1  and  2  is  50%  and 
921%,  respectively.  By  processing  three  data  sets,  the  RSS  analysis  yields  solutions 
for  the  two  parameters  with  an  error  of  0.2%  and  3%,  respectively.  With  the  dras¬ 
tic  improvement,  a  sample  estimation  case  (Case  IA)  is  attempted  using  the  most 
accurate  RSS  result.  This  estimation  case  converged  and  is  therefore,  considered 
a  success.  However,  with  no  noise  and  with  a  near-perfect  guess  for  the  nonlinear 

dynamics,  one  would  have  expected  a  nearly  perfect  estimate  with  a  distribution  of 
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residuals  closer  to  zero.  This  was  not  the  outcome.  In  fact,  had  <7instr  been  set  to 
a  value  closer  to  machine  zero,  the  standard  deviation  criterion  would  have  failed. 
This  was  discovered  by  allowing  the  estimator  to  continue  its  iterations  past  con¬ 
vergence  where  aresid  remained  practically  constant  at  5.77835  x  10-5  from  iteration 
56  to  200.  Several  other  attempts  were  made  using  different  values  for  k\  with  no 
significant  change  in  the  outcome.  It  appears  the  estimator  reached  a  limit  of  some 
type  preventing  any  improvement  of  the  estimated  solution.  This  limit  seen  in  the 
standard  deviation  remaining  nearly  constant  for  144  iterations  is  most  likely  due 
to  the  estimator  finding  a  local  minimum  on  a  nearly  level  surface  closely  shared 
with  the  desired  solution.  By  adjusting  the  initial  conditions  of  the  relative  position 
states  at  epoch  to  values  closer  to  the  truth  quickly  resolved  this  problem.  Case  IA 
simply  shows  the  sensitivity  of  the  estimator  and  the  outcome  based  solely  on  the  a 
priori  values  used  to  start  the  estimation  process. 

4-5.2  Case  II.  Case  II  (12-State),  which  used  a  significantly  larger  data  set 
(6053  measurements),  proved  to  be  a  failed  estimation  attempt  without  convergence. 
As  mentioned  in  the  results,  <Jrestlj  remained  at  14.042  meters  from  iteration  28  to  200. 
The  result  of  this  can  be  seen  in  the  radical  rotational  dynamics  behavior.  From 
Figure  4.13,  it  appears  that  the  singularity  in  the  Euler  angles  is  contributing  to 
the  uncharacteristic  rotational  dynamics  based  on  the  incorrect  reference  solution  at 
epoch.  The  existence  of  multiple  solutions  is  a  valid  concern  here  as  it  was  with  Case 
IA,  where  an  accuracy  limit  was  reached  by  the  estimator.  The  noisy  data  simply 
introduces  more  difficulty  for  the  estimator.  In  this  case,  the  incorrect  reference 
solution  is  caught  by  the  convergence  criteria;  specifically,  the  distribution  of  the 
residuals.  The  9-State  portion  of  Case  II  illustrated  the  multi-solution  problem 
quite  well.  All  three  trials  quickly  converged  to  a  solution  in  8  to  9  iterations  in 
under  120  seconds  with  aresid  ~  crirastr. .  However,  the  difference  is  best  seen  in  the 
nonlinear  rotational  dynamics.  The  rotational  equations  of  motion  are  nonlinear 
and  are  extremely  sensitive  to  the  initial  conditions.  This  is  illustrated  in  the  figures 
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propagating  the  dynamics.  The  estimation  would  lead  one  to  believe  that  the  RSO  is 
in  a  single-axis  spin.  The  error  in  the  estimated  MOI  values  are  clearly  related  to  the 
error  of  the  a  priori  initial  estimates.  This  highlights  the  known  fact  in  estimation, 
that  the  better  the  initial  guess,  the  better  the  estimate. 

4-5.3  Case  III.  By  introducing  a  second  data  batch  to  the  estimation 
process  for  Case  III  (12-State),  the  estimation  results  radically  improved.  Using 
the  same  number  of  data  points  per  noisy  data  set  as  in  Case  II,  the  results  of 
this  case  far  exceed  that  of  Case  1A  and  certainly  Case  II  (12-State)  in  terms  of 
accuracy.  In  100  seconds  and  13  iterations  convergence  was  reached  resulting  in  an 
estimated  solution  with  ar esid  within  0.4%  of  <Jmstr  for  the  first  data  set  and  1.6%  for 
the  second  data  set.  In  the  9-State  estimation  runs,  trials  one  and  two  successfully 
converged  in  18  and  21  iterations,  respectively  and  required  213  seconds  and  268 
seconds  of  processing  time,  respectively.  The  rotational  dynamics  slightly  improved 
while  the  moments  of  inertia  suffered.  Trial  3  failed  the  standard  deviation  criterion 
and  therefore,  did  not  converge  on  a  solution  in  under  200  iterations  (48  minutes). 
In  fact,  the  Trial  3  run  reached  a  constant  aresid  =  0.32618  meters  (for  the  n  =  2 
data  set)  at  iteration  76  until  program  termination  at  iteration  200.  With  two  data 
sets  used  in  this  estimation  attempt  and  after  convergence  was  reached  for  this  trial 
(using  the  same  a  prior  information)  in  Case  II,  it  is  unlikely  that  solution  ambiguity 
is  the  primary  fault.  The  other  explanation  for  this  failed  attempt  may  be  that  the 
moments  of  inertia  are  not  fully  observable.  This  possible  cause  is  first  apparent  in 
this  case. 

4-5-4  Case  IV.  In  Appendix  C,  the  results  of  a  fourth  case  are  provided. 
This  extra  case  uses  the  same  initial  conditions  and  filter  settings  as  the  previous  12- 
State  cases;  however,  three  batches  of  data  are  processed  together.  The  results  of  this 
case  are  closer  to  the  true  values  than  Case  III  (12-State),  emphasizing  the  improved 
accuracy  using  multiple  data  sets  in  the  estimation  process  for  the  12-State  vector 
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version.  However,  the  outcome  was  not  successful  for  the  9-State  estimation  cases. 
Using  the  same  settings  as  in  Case  II  and  III  (9-State),  and  making  use  of  a  third 
data  set,  the  estimator  did  not  achieve  convergence  within  the  specified  criteria  and 
all  three  trials  failed.  Details  of  these  runs  are  found  in  Appendix  C.  All  three  trials 
failed  as  a  result  of  the  residual  distribution  of  the  second  data  set  (associated  with 
n  =  2).  With  (i-resid  ~  0.327,  similar  to  Case  III  (9-State),  the  convergence  criterion 
was  not  satisfied  for  any  of  the  three  trials  in  200  iterations.  It  is  interesting  to  note 
that  convergence  failure  was  attributed  to  data  set  two  and  that  aresid  for  all  three 
trials  were  nearly  equal.  This  is  the  first  case  where  the  introduction  of  an  additional 
data  set  actually  resulted  in  poor  performance  and  estimation  failure  for  all  trials. 
With  three  data  sets,  it  is  difficult  to  associate  the  failure  mechanism  completely  to 
solution  ambiguity.  In  fact,  with  this  case  it  is  more  apparent  that  the  moments  of 
inertia  may  simply  not  be  observable. 

4-5.5  Comparison  of  12-State  Cases.  Table  4.24  compares  the  12-State 
estimate  errors  in  terms  of  the  number  of  data  sets  processed.  As  mentioned  in  the 
results,  the  12-State  estimate  using  one  data  set  took  200  iterations  until  program 
termination.  The  12-State  estimation  runs  using  two  and  three  data  sets  converged 
in  13  iterations  (approximately  130  seconds).  The  table  illustrates  the  improvement 
in  estimation  accuracy  as  the  number  of  data  sets  used  is  increased  and  solution 
ambiguity  is  decreased. 

4-5.6  Comparison  of  9-State  Cases.  A  similar  comparison  can  be  made 
using  Tables  4.17,  4.23,  and  the  results  presented  in  Appendix  C.  This  section  will 
compare  each  9-State  estimation  trial  separately.  For  each  trail,  a  comparison  of  the 
number  of  data  sets  used  and  the  resulting  error  will  be  made.  For  Trial  1,  the  a 
priori  estimates  for  the  MOI  have  an  initial  error  of  10%  from  the  true  values.  The 
average  MOI  errors  (the  average  of  components  A,  B,  and  C  together)  for  runs  made 
using  1,  2,  and  3  data  sets  are  11%,  26%,  and  37%,  respectively.  In  this  case,  the 
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Table  4.24  12-State  Estimation  Error  Summary 


Estimate  Error 

State 

1  Data  Set 

2  Data  Sets 

3  Data  Sets 

[Case  II] 

[Case  III] 

[Case  IV] 

5x 

[m] 

35.9% 

0.02% 

0% 

8y 

[m\ 

204.9% 

0.08% 

0% 

6z 

[m] 

52.5% 

0.1% 

0.0002% 

8x 

[m/s] 

122% 

0% 

0% 

Sy 

[m/s] 

309% 

0% 

0% 

5z 

[m/s] 

59% 

0% 

0% 

CJl 

[rad/ s] 

0.0148 

0.0000 

0.0000 

00  2 

[rad/ s] 

0.0434 

0.0000 

0.0000 

Ll>3 

[rad/ s] 

0.0392 

0.0000 

0.0000 

0i 

[rad] 

4.7289 

0.0006 

0.0001 

62 

[rad] 

0.9465 

0.0004 

0.0001 

03 

[rad] 

8.5974 

0.0008 

0.0018 

MOI  error  using  a  single  data  set  is  close  to  the  a  priori  error  of  10%.  Otherwise, 
increasing  the  number  of  data  sets  processed  increased  the  moments  of  inertia  error 
overall  while  the  error  for  the  rotational  states  slightly  decreased. 

For  Trial  2,  the  initial  a  priori  estimates  for  the  MOI  have  an  error  of  20%. 
In  this  case,  the  average  MOI  error  for  estimation  runs  made  with  1,  2,  and  3  data 
sets  are  20%,  17%,  and  18%,  respectively.  Increasing  the  number  of  data  batches 
processed  did  not  improve  the  moments  of  inertia  estimates  yet  it  slightly  improved 
the  rotational  dynamics  estimates.  The  average  MOI  error  is  relatively  close  to  the 
a  priori  MOI  error. 

Trial  3  uses  initial  values  for  the  moments  of  inertia  which  are  30%  in  error 
from  the  true  values.  In  this  Trial  the  two  runs  using  multiple  data  sets  (from  Cases 
III  and  IV)  failed  to  converge.  The  average  MOI  errors  for  estimation  runs  involving 
1,  2,  and  3  data  sets  are  30%,  34%,  and  34%.  Once  again,  the  average  error  is 
comparable  to  the  initial  error  used  to  begin  the  estimation  process. 
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4-6  Summary  of  Results 

From  the  data  presented  here,  it  can  be  seen  that  multiple  solutions  that  satisfy 
the  convergence  criteria  do  in  fact  exist  for  this  scenario.  The  solution  uncertainty 
can  be  significantly  reduced  while  improving  solution  accuracy  for  the  estimation 
of  the  12-State  vector  by  processing  multiple  data  sets  (as  shown  in  the  12-State 
estimation  cases).  It  can  be  said  that  accurate  estimation  of  the  relative  motion 
and  rotational  states  proved  to  be  possible.  The  direct  estimation  of  the  moments 
of  inertia  however,  proved  to  be  more  of  a  challenge  with  relatively  lower  accuracy 
than  the  other  estimated  states  if  the  convergence  criteria  was  even  achieved.  This 
may  be  partially  due  to  the  existence  of  multiple  solutions  but  more  importantly,  the 
lack  of  observability  of  the  moments  of  inertia  is  the  likely  cause.  If  the  moments  of 
inertia  are  not  fully  observable  based  on  range  measurements,  the  estimator  would 
then  not  be  fully  capable  of  accurately  estimating  the  MOI  (based  on  the  equations 
of  motion)  yielding  poor  results  for  the  9-State  solutions.  In  such  a  case,  processing 
multiple  data  sets  affects  only  solution  ambiguity  but  does  not  improve  observability. 
This  explanation  best  fits  the  results  seen  by  this  research  for  the  9-State  vector 
estimation  runs.  Slight  improvements  in  the  rotational  dynamics  are  seen  which 
can  be  attributed  to  the  reduction  of  multiple  solutions  using  the  additional  data 
sets.  Unlike  the  improving  rotational  state  estimates,  the  additional  data  sets  do  not 
improve  the  MOI  estimates.  This  conclusion  will  be  discussed  further  in  the  next 
Chapter.  Note  that  for  all  of  the  MOI  estimates  (successful  and  unsuccessful),  the 
components  are  all  positive  and  they  all  satisfy  the  true  MOI  relation  of  A  >  B  >  C . 
Therefore,  although  the  individual  components  of  the  moment-of- inertia  matrix  are 
not  well  estimated,  the  proper  relationship  associated  with  the  major,  minor,  and 
intermediate  axes  of  inertia  of  the  simulated  Resident  Space  Object  are  maintained. 
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V.  Conclusions  &  Recommendations 

5. 1  Summary 

This  thesis  demonstrates  the  basic  ability  to  accurately  estimate  the  relative 
motion  and  rotational  dynamics  of  a  simulated  Resident  Space  Object  using  a  non¬ 
linear  least  squares  estimation  filter  based  on  range  observations  collected  during  a 
proximity  mission.  However,  the  feasibility  of  accurately  estimating  the  moments  of 
inertia  properties  of  the  RSO  is  questionable  and  proved  to  be  challenging.  In  fact, 
this  research  concludes  that  the  moments  of  inertia  of  the  RSO  may  not  be  fully 
observable  using  range  observations  alone.  This  thesis  indicates  that  processing 
multiple  data  batches  significantly  reduces  solution  ambiguity  and  improves  estima¬ 
tion  accuracy,  as  seen  with  the  12-State  vector  estimation  cases.  Conclusions  from 
these  estimated  solutions  can  be  drawn  to  further  improve  the  understanding  of  the 
dynamical  and  physical  properties  of  the  RSO. 

5.2  Conclusions 

This  research  indicates  that  implementing  a  nonlinear  least  squares  estimator 
can  successfully  result  in  an  estimated  solution  for  the  relative  motion  dynamics, 
rotational  dynamics,  and  moments  of  inertia  of  an  RSO  with  certain  limitations. 
The  results  from  Case  1A  and  Case  II  (12-State)  illustrate  the  critical  role  the  initial 
a  priori  estimates  play  in  the  overall  success  or  failure  of  the  estimation  attempt. 
The  existence  of  multiple  estimation  solutions  may  prove  to  be  quite  problematic 
in  the  case  of  an  uncooperative  RSO.  One  must  discern  between  multiple  solutions 
that  are  all  mathematically  valid  and  the  single  solution  representative  of  the  true 
states.  This  is  challenging  when  the  luxury  of  knowing  the  truth  is  not  always 
realistic.  However  taking  multiple  range  measurements  of  tracked  points  on  the  RSO 
body  provides  multiple  data  sets.  Each  data  set  significantly  reduces  the  solution 
ambiguity  as  shown  in  the  RSS  surface  plots  and  demonstrated  in  Case  III  and  IV 
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(12-State).  This  can  further  enhance  one’s  estimation  attempt  and  increase  solution 
believability  and  accuracy. 

In  this  research  the  limitation  of  estimating  the  dynamical  properties  of  the 
RSO  is  limited  to  the  deterministic  two  body  motion  developed  in  Chapter  2.  This 
is  well  suited  for  estimating  the  dynamics  where  the  RSO  is  a  small  heavenly  body 
(e.g.  an  asteroid).  However,  when  the  RSO  is  a  manmade  spacecraft,  this  quickly 
becomes  a  complex  problem.  RSO  motion  not  fully  described  by  classic  two-body 
dynamics  would  not  be  suited  for  the  estimator  used  in  this  research.  Thruster  brings 
for  maneuvering  or  station-keeping,  energy  dissipation,  and  external  perturbations 
are  all  realistic  contributors  to  the  overall  RSO  dynamics,  which  are  not  captured 
by  this  ideal  model.  From  the  model  used  in  this  research,  and  from  the  simulation, 
one  can  say  that  the  simulated  RSO  is  rotating  primarily  about  two  axes  at  certain 
estimated  angular  velocities  (uq  and  uj3).  At  any  instant  in  time,  one  could  predict 
the  orientation  of  the  RSO  body.  Also,  with  some  knowledge  of  the  moments  of 
inertia,  the  geometry,  and  the  estimated  angular  velocities,  one  can  conclude  that 
the  RSO  is  spinning  primarily  about  its  major  axis  of  inertia  and  is  therefore,  in 
a  relatively  stable  state.  This  study  concludes  that  the  relative  position,  relative 
velocity,  angular  velocity,  and  orientation  angle  states  can  be  accurately  estimated 
using  the  12-State  vector  version  by  processing  multiple  range  data  sets. 

Perhaps  a  more  challenging  task  is  to  gain  knowledge  of  the  physical  properties 
of  a  RSO.  This  is  true  for  both  heavenly  bodies  as  well  as  uncooperative  Targets.  At 
a  very  basic  level,  this  research  shows  that  some  knowledge  of  the  physical  properties 
must  be  known  or  estimated  beforehand.  The  RSO  center-of-mass  location  is  of  key 
importance  to  the  following:  the  CW  equations  (for  relative  motion  dynamics),  the 
coordinate  frame  transformations,  the  observational  points  measured  from  the  CM, 
and  finally,  for  the  moment-of-inertia  components  themselves.  The  9-State  cases  in 
this  thesis  demonstrate  that  the  direct  estimation  of  the  moments  of  inertia  using 
range  measurements  alone  is  problematic.  This  research  concludes  that  the  moments 
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of  inertia  may  not  be  fully  observable  using  the  current  dynamics  model  and  range 
observations.  In  this  case,  the  use  of  multiple  data  sets  simply  improves  solution 
ambiguity  (as  seen  in  the  rotational  states)  but  does  not  affect  observability  (as  seen 
in  the  error  of  the  MOI  estimates).  The  quality  of  the  MOI  estimated  results  proved 
to  be  directly  related  to  the  quality  of  the  a  priori  initial  guess.  Therefore,  a  more 
detailed  study  is  required  to  examine  the  observability  conditions  and  expressions 
needed  to  accurately  and  directly  estimate  the  moments  of  inertia. 

5. 3  Contributions 

In  the  held  of  spacecraft  proximity  operations,  this  research  illustrates  the 
application  and  performance  of  a  nonlinear  least  squares  estimation  filter  for  relative 
motion  and  rotational  dynamics  using  range  observations  made  by  a  spacecraft  in  the 
vicinity  of  a  RSO.  Perhaps  of  more  significance,  this  study  exposes  the  limitations 
of  the  implemented  dynamics  model  and  estimator  in  regard  to  the  observability  of 
the  moments  of  inertia  of  the  Resident  Space  Object. 

5.4  Recommendations 

There  are  several  areas  in  which  additional  research  would  be  worthwhile  to 
pursue  with  this  baseline  ability  established.  Continuing  with  this  nonlinear  least 
squares  estimation  approach,  it  would  be  of  benefit  to  employ  some  realism  to  both 
the  truth  model  and  the  data  sets.  This  would  include  incorporating  J2  perturba¬ 
tions,  solar  radiation  pressure,  and  perhaps,  atmospheric  drag.  To  make  the  truth 
model  robust  and  versatile  it  would  be  beneficial  to  use  quaternions  to  describe  the 
Target’s  orientation  and  avoid  singularity  limitations.  This  research  used  data  sets 
which  have  data  points  at  constant  time  intervals.  It  would  be  of  interest  to  study 
the  effects  of  using  data  sets  with  measurements  not  taken  entirely  at  constant  time 
intervals.  This  would  simulate  an  intermittent  data  acquisition  stream  as  well  as 
blackout  periods  when  the  tracked  point  is  not  visible  to  the  sensor  due  to  physi- 
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cal  obstruction.  This  leads  to  additional  research  in  the  amount  of  data  needed  for 
successfully  estimating  the  true  states.  LIDAR  sensors  which  are  now  being  used 
frequently  in  rendezvous  and  proximity  operations  can  provide  data  in  a  range  and 
two  angle  set  or  in  a  range  vector  set  [12].  Introducing  directional  information  into 
the  data  would  certainly  refine  the  estimate  and  also  reduce  solution  ambiguity.  Ad¬ 
ditionally,  researching  how  the  selected  locations  of  the  observational  points  on  the 
RSO  affect  the  estimator’s  results  would  prove  beneficial.  Of  significant  importance, 
research  should  be  conducted  into  the  observability  of  the  moments  of  inertia. 

It  would  of  interest  to  incorporate  an  extended  Kalman  filter  into  this  estima¬ 
tion  process.  This  would  simply  involve  a  rearrangement  and  modification  of  the 
existing  matrices  in  the  estimator  code.  Using  at  least  two  data  sets  for  some  batch 
size,  a  least  squares  estimation  is  performed  and  at  some  point  handed  over  to  the 
Kalman  filter  to  continue  the  estimation  of  the  state  variables.  This  would  be  very 
useful  in  the  12-State  variable  case. 
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Appendix  A.  Linearized  Observation 

This  appendix  contains  the  calculations  used  to  develop  the  partial  derivatives  which 
are  used  to  calculate  the  H  matrix.  Beginning  with  the  derived  range  observation 
function 

G  =  rn=  yj (pXn  -  5x )2  +  {pVn  -  Sy )2  +  (pZn  -  ~Sz)2  (A.l) 

where  the  subscript  n  denotes  observation  point  number  n.  Recall  the  definition 

<9G 

HdL)  =  (A. 2) 


Note  that  the  range  equation  above  is  a  function  of  Sx,  Sy.  Sz,  pXn ,  pVn,  and  pZn. 
In  terms  of  the  dynamics  state  variables,  the  range  equation  is  only  a  function  of 
8x,  Sy,  Sz,  6\,  62,  and  #3.  For  simplicity,  all  other  state  variables  which  the  range 
equation  is  not  a  function  of  are  ignored  till  the  end.  This  results  in 


H  = 


dG 

dX 


dG 

dSx 


dG 

dSy 


dG 

dSz 


dG 

de  1 


dG 

d02 


dG 

de3 


(A.3) 


Recalling  that  pXn ,  pVn ,  and  pZn  are  functions  of  the  Euler  angles  due  to  the  rotation 
of  the  Target  body,  the  H  matrix  can  simply  be  found  using  the  chain  rule  in  the 
following  fashion 


H 


dG  dG  dG  dG 

dSx  dSy  dSz  dpXn 


1 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

dG 

dG 

0 

0 

1 

0 

0 

0 

dPyn 

dPzn. 

0 

0 

0 

9pxn 

de  1 

9Pxn 

de2 

9px„ 

d$3 

0 

0 

0 

dPVn 

001 

dPyn 

do2 

dPyn 

dOs 

0 

0 

0 

9Pzn 

de  1 

dPzn 

d02 

1 

M  00 

(A.4) 
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To  compute  the  partial  derivatives  recall  that 


llr'll  =  ||[R-'6]Et.^  -  <5,0*11  (A. 5} 

where  the  position  vector  in  the  orbital  frame  of  the  observation  point  is  defined  as 

(A. 6) 

define  the  position  vector  in  $e  as 


Rf:  =  [Re6lR^ 


Rn 


define  the  position  vector  in  as 


(A.7) 


(A.8) 


Then  substituting  the  vector  definitions  into  Equation  A. 6  yields 

C2C3  S1S2C3  —  C1S3  C1S2C3  +  S1S3 

c2s3  S1S2S3  +  cic3  C1S2S3  -  S1C3 

—  S2  S1C2  C1C2 

carrying  out  the  matrix  multiplication  results  in  the  following  set  of  equations 


(A.9) 


|  Pxn 

^1„C2C3  +  62n(SlS2C3  -  Ci S3)  +  63„(ClS2C3  +  S1S3) 

|  Pyn 

>  = 

frlnC2S3  +  &2n(siS2S3  +  C1C3)  +  63rl(ciS2S3  -  S1C3) 

U*, 

-fri„s2  +  b2n  S!C2  +  &3„CiC2 

(A.10) 
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Taking  the  partial  derivatives  of  each  equation  with  respect  to  6\,  d2,  03  results  in 


dPxn 
00! 
O 'Pxn 
00  2 
9Pxn 

00 3 

dPyn 
O0X 
dPvn 
00 2 
dPvn 
00  3 
OPzn 
00! 
9Pzn 
002 


&2ns2c3ci  +  62„s3Si  -  63ns2c3Si  +  63nCiS3  (A. 11) 

^inc3s2  +  62nSic2c3  +  b3n  Cic2c3  (A. 12) 

-&i„c2s3  -  b2n sis2s3  -  62„cic3  -  63„cis2s3  +  63„sic3  (A. 13) 

^2ncis2s3  -  62nsic3  -  b3n sis2s3  -  63ncic3  (A. 14) 

-&i„s2s3  +  62„sic2s3  +  b3n  cic2s3  (A. 15) 

^i„c2c3  +  62n sis2c3  -  62ncis3  +  b3n cis2c3  +  63nsis3  (A. 16) 

&2„cic2  -  63„sic2  (A. 17) 

-^inc2  -  &2nsis2  -  53ncis2  (A. 18) 


it  =  0  (A'19) 

Now  taking  the  partial  derivatives  of  the  range  function,  Equation  A.l,  with  respect 
to  Sx,  Sy,  and  Sz  results  in 


OG 

Sx  -  pXn 

OSx 

V ( Pxn  -  Sx)2  +  ( pyn  -  5y)2  +  (pZn  -  Sz) 

OG 

Sy  ~  Pyn 

OSy 

V ( Pxn  -  Sx)2  +  (pyn  -  5y)2  +  0 pZn  ~  Sz) 

OG 

Sz  -  pZn 

05  z 

\J ( Pxn  ~  Sx)2  +  (pyn  -  5y)2  +  ( pZn  -  Sz) 

(A. 20) 
(A.21) 
(A. 22) 


These  12  partial  derivatives  are  coded  into  a  MATLAB  function  file  in  the  form  of 
Equation  A. 4.  The  computer  carries  out  the  matrix  multiplication,  and  by  the  chain 
rule,  solves  for  the  H  matrix.  Depending  on  which  version  of  the  state  matrix  is 
used,  the  program  forms  the  proper  version  of  the  H  matrix. 
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For  the  12-state  vector  version  defined  as 


X  = 


1  T 


Sx  Sy  Sz  Si:  Sy  Sz  uj\  uo 2  UJ3  6\  62  O3 


the  linearized  observation  function  is  the  solution  to 


H 


<9G  dG  qg  n  n  n  n  n  n  aG  9G  9G 

dSx  d8y  d8z  U  U  U  U  U  U  ddi  dd2  d03 


For  the  9-state  vector  version  defined  as 


X  = 


l  T 


uj\  ui2  ^3  9\  62  O3  A  B  C 


the  linearized  observation  function  is  the  solution  to 


H  = 


000^^^000 

u  u  u  96»1  d02  de3  u  u  u 
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(A. 23) 


(A. 24) 


Appendix  B.  Linearized  Dynamics 

Recall  an  analytical  solution  to  the  Clohessy-Wiltshire  equations  exists.  This  so¬ 
lution  was  presented  in  the  form  of  a  state  transition  matrix,  <&cw,  in  Chapter  2. 
Therefore,  the  state  transition  matrix  for  the  rotational  portion  of  the  dynamics  is 
in  need.  As  defined  in  Chapter  2,  the  rotational  state  vector  is 


Xrot(£)  — 


OJl  ii>2  OJ3  61  62  O3 


(Bl) 


The  time-derivative  of  the  rotational  state  vector  is  then  given  below.  The  vector 
function,  f ,  is  then  associated  with  each  state  equation  of  motion. 
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A 

(B-2) 


The  A  matrix  is  defined  as  the  matrix  containing  the  partial  derivatives  of  the 
equations  of  motion  with  respect  to  the  state  variables. 


dfi  dfi  dfi  dfi  dfi  dfi 
dw!  duo2  duo3  dO 1  d02  d03 


d/2  d/2  d/2  d/2  d/2  d/2 
duo  1  duo 2  duo 3  dOi  d02  d03 


d/3  d/3  d/3  d/3  d/3  d/3 

^-11 

q3x3 

<9o;i  du;2  do;3  <9#i  <902  d#3 

d/4  d/4  d/4  d/4  d/4  d/4 

A21 

A22 

<9o;i  <9a;2  du;3  <9#i  d#2  d$3 

d/5  d/5  d/5  d/5  d/5  d/5 
duo  1  duo 2  duo 3  dO  1  d02  d03 


d/6  d/6  d/6  d/6  d/6  d/6 

<9o;i  du;3  <9#i  <902  d#3 
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The  A  matrix  is  partitioned  in  for  clarity.  The  partial  derivatives  comprising  each 
sub-matrix  are  given  as 


0 

B-C.  . 

A  ^3 

B-C 

A 

C-A,  , 

B 

0 

C-A 
B  ' 

a-b ,  . 

C 

A-B,  , 

C  UJl 

0 

1  Sit2  Cit2 
0  ci  —si 

0  Si sec2  Ci sec2 

OJ  2CA2  —  UT^Si 
— CU2S1  —  W3C1 


u^Si  sec^  +W3C1  sec^  0 


LJ2C1  sec2  —  W3S1  sec2  CU2S1  sec2 12  +  W3C1  sec2 12  0 
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where  s,,:  =  sin  9X ,  cx  =  cos  9X ,  tx  =  tan  0X .  secx  =  sec  9X ,  and  sec2  =  sec2dx.  These 
partial  derivatives  are  coded  into  a  MATLAB  function  file  for  the  A  matrix. 

The  state  transition  matrix  for  the  rotational  portion  is  defined  as 


011 

012 

013 

014 

015 

016 

021 

022 

/7)00 

023 

(T)o  0 
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(1)0  A 
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and  the  time-derivative  of  the  state  transition  matrix  is  simply  defined  as 


011 

012 
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014 
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The  matrices  above  are  coded  into  the  estimation  program  and  using  a  Runga-Kutta 
algorithm,  the  numerical  solution  to 


&rot(t,t0)  =  A(t)$  rot  (Mo) 


(B.9) 


is  solved  for  resulting  in  $>rot(t,to).  Incorporating  the  CW  portion  results  in  the 
12  x  12  state  transition  matrix  for  the  12-state  vector  version  used  in  the  estimation 
process 


d>6x6 

CW 

q6x6 

q6x6 

<E>6x6 

^ rot 

(B.  10) 


In  a  similar  fashion,  the  9-state  vector  version  of  the  A  matrix  and  the  state 
transition  matrix  is  found.  This  section  should  be  treated  separate  from  the  12- 
state  version  above.  Re-defining  similar  variables,  omitting  the  CW  portion  of  the 
state  vector,  and  including  the  moment-of-inertia  variables,  the  9-state  version  is 
formed.  Beginning  with  the  9-state  vector  defined  as 


X  = 


1  T 


(jj\  UO2  ^3  9\  62  O3  A  B  C 
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the  time-derivative  of  the  state  variables  are  taken  resulting  in 
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The  9-state  version  of  the  A  matrix  is  defined  as 


>fi  dfi  dfi  dfi  dfi  dfi  dfi  dfi  dfx 
duoi  duj2  8003  d$i  862  863  8 A  8B  8C 

8f 2  d/2  d/2  d/2  d/2  d/2  d/2  d/2  df2 
duoi  8uj2  duo 3  86i  802  803  8 A  8B  8C 

8f 3  8f3  8 f3  df3  8 f3  df3  <9/3  8f3  df3 
duoi  8uj2  duo 3  80 1  802  80 3  8 A  8B  8C 

8/4  8f 4  d/4  8f 4  d/4  d/4  8f 4  8f 4  8f 4 

duji  duo 2  8003  80 1  802  80 3  8 A  8B  8C 

a  6x6 
^ rot 

^3x3 

A  = 

8f^8f^8f^8h8h8h8h8h8h 

duoi  duo2  8003  80 1  802  80 3  8 A  8B  8C 

= 

q3x3 

df6  df6  df6  df6  df6  df6  df6  df6  df6 
duoi  duo2  duo 3  80i  80 2  80 3  8 A  8B  8C 

q3x9 

8fr  dfr  8fj_  dfr  dfr  dfr  dfr  dfr  dfr 
duoi  8002  8003  80 1  80 2  80 3  8 A  8B  8C 

d/s  d/s  df8  df8  d/8  df8  df8  df8  d/8 
duoi  duo2  8003  80 1  80 2  80 3  8 A  8B  8C 

df9  d/9  d/9  d/9  d/9  d/9  d/9  df9  df9 

_8uoi  8002  8003  80i  80 2  80 3  8 A  8B  8C  _ 

The  sub-matrix  A®*t6  is  further  partitioned  into 


A  6x6 

^ rot 


q3x3 

A21 

A22 

B-4 

(B.  11) 


(B.  12) 


(B.  13) 


where  sub-matrices  An,  A2 1,  A22,  are  equal  to  Equations  B.4,  B.5,  and  B.6,  respec¬ 
tively.  The  solution  to  sub-matrix  B  is 


B6XS  = 


(B-C) 

'  ^2  UJ2Ll)3 

^2<^3 

(C-A),  ,  ,  , 

^1^3 

—  -£jCUlUJ2  — 

~^UJ2(jJ3 

jjLOlUJS 
( A-B) 


(B.14) 


This  completes  the  solution  to  the  9x9  matrix  A  .  The  9x9  state  transition  matrix 
is  then  defined  as 


011  012  013 
021  022  023 
031  032  033 

041  042  043 

^  =  051  052  053 

061  062  063 
071  072  073 
081  082  083 
091  092  093 


014  015  016  017  018 

024  025  026  027  028 

034  035  036  037  038 

044  045  046  047  048 

054  055  056  057  058 

064  065  066  067  068 

074  075  076  077  078 

084  085  086  087  088 

094  095  096  097  098 


019 

029 

039 

049 

059  (B.15) 

069 

079 

089 

099 
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and  the  time-derivative  of  the  state  transition  matrix  is  defined  as 

011  012  013  014  015  016  017  018  019 

021  022  023  024  025  026  027  028  029 

031  032  033  034  035  036  037  038  039 

041  042  043  044  045  046  047  048  049 

051  052  053  054  055  056  057  058  059  (B.16) 

061  062  063  064  065  066  067  068  069 

071  072  073  074  075  076  077  078  079 

081  082  083  084  085  086  087  088  089 

091  092  093  094  095  096  097  098  099 

The  A,  <l>,  and  4?  matrices  above  are  coded  into  the  estimation  program  and  using 
a  Runga-Kutta  numerical  integration  algorithm,  the  numerical  solution  to 

t0)  =  A(t)®(t,  t0)  (B.17) 

is  solved  for.  This  results  in  the  numerical  solution  of  the  9x9  state  transition 
matrix  d>(t,f0)  for  the  9- variable  state  vector. 
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Appendix  C.  Additional  Data 


Truth  Model  Data 

The  Truth  Model  program  settings  and  constants  are  shown  in  Figure  C.17 


CONSTANTS  £  USER  PARAHETERS 


mu,  Earth  =  393600.44  [kmA3/secA2]  Radius  of  Earth,  Re  =  6373.14  [km] 

Mean  Orb  Motion,  Target  =  1.0331e-003  [rad/s]  Orbital  Period  =  6052  [s] 

Target  Altitude  =  300.00  [km] 


Tgt  Orbit  Incl 

Assumed  Target 

MO  I :  A 

=  23.50  [deg] 

=  45.OO00 

Tgt  Orbit 

B  =  30.0000 

RAAN  =0.00  [deg] 

C  =  15.0000 

INITIAL  STATE 

CONDITIONS 

POSITION  [m] 

VELOCITY  [m/s] 

ANGULAR  VEL  [rad/s] 

EULER  ANGLE  [rad] 

dx=  5O.0O0O 

dy=  5O.0O0O 

dE=  50.0000 

dx  do t= 

dy  do t= 

ds  dot= 

0.0100 

0.0100 

0.0100 

omegal=  O.OIOOOOOO 

omega2=  O.OOOOOOOO 

omega3=  0.00100000 

thl=  0.00000000 

th2  =  0.00000000 

th3=  0.00000000 

#  of  ORBITAL  PERIODS 

START  TIME  [s]  END  TIME  [s] 

Int  Step-Sise  [s] 

1.0000 

0.0000 

6052 .0000 

1.0000 

Figure  C.l 

Truth  Model 

Case  I  Additional  Data 

Additional  data  from  Case  1A  are  provided  in  the  figures  below. 


Case  II  Additional  Data 
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Relative  Spacecraft  Position  Error 


Figure  C.2  Case  I:  True  Error  of  Relative  Position  Estimates 
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Figure  C.3  Case  I:  Estimated  Relative  Velocity  States 


Figure  C.4  Case  I:  True  Error  of  Relative  Velocity  Estimates 
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Figure  C.5 


Case  I:  Estimated  Angular  Velocity  States 


Angular  Velocity  Error 
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Figure  C.6 


Case  I:  True  Error  of  Angular  Velocity  Estimates 


Euler  Angle  Orientation  Error 


Figure  C.7  Case  I:  True  Error  of  Euler  Orientation  Angle  Estimates 
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Figure  C.8  Case  II:  True  Error  of  Relative  Position  Estimates(12-State) 


Figure  C.9  Case  II:  True  Error  of  Relative  Velocity  Estimates  (12-State) 
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Figure  C.10  Case  II:  True  Error  of  Angular  Velocity  Estimates  (12-State) 
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Figure  C.ll  Case  II:  True  Error  of  Euler  Orientation  Angle  Estimates  (12-State) 


Actual  Range  Measurements  vs.  Predicted  Range  Measurements 
Target  Observation  Point  1 


Figure  C.12  Case  II:  Range  Profile  Close-Up  View  (9-State) 
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Error  8  z  -vel  [m/s]  Error  8  y  -vel  [m/s]  Error  8  x  -vel  [m/s] 


Case  III  Additional  Data 


Figure  C.13  Case  III:  True  Error  of  Relative  Position  (12-State) 


Figure  C.14  Case  III:  True  Error  of  Relative  Velocity  (12-State) 
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Figure  C.15  Case  III:  True  Error  of  Angular  Velocities  (12-State) 
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Figure  C.16  Case  III:  True  Error  of  Euler  Orientation  Angles  (12-State) 
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Case  IV  Data 


RANGE  DATA  TO  PROCESS 


Number  of  Observations  =  6053  Obs  Time-Step  =  1.00  [3] 

Start  Epoch  Obs  lime  =  0.00  Is]  End  Obs  Time  =  6052.00  [s] 


FILTER  PERFORMANCE 


Convergence  Status :  DID  CONVERGE 

Filter  Terminated  on  Iteration  Number:  13 

Standard  Deviation  of  Data  Residuals  from  Final  Iteration:  0.2479376954 

Standard  Deviation  of  Data  Residuals  from  Final  Iteration:  0.2496057646 

Standard  Deviation  of  Data  Residuals  from  Final  Iteration:  0.2522780251 

CPU  Processing  Time  =  129.26714852  [3] 

*******************  it  *********  ******  4***4  *****  ******  *****  ****************  ********** 

N-LSQ  ESTIMATOR  RESULTS 

*****************************************************************4**************** 

POSIT  I  OH  [id]  VELOCITY  [m/s]  ANGULAR  VEL  [rad/s]  EULER  AHGLE  [rad] 


dx=  50.00097336 
dy=  50.00165926 
dz=  49.99363954 


d>:dot=  0.0100 02  8  6 
dydot=  0.00999058 
dtdot=  0.00996686 


omgl=  0.00999977 
omg2=  -0.00001785 
omg3=  0.00100841 


thl=  -0.00007261 
th2=  0.00009480 
th3=  0.00176362 


Version:  12-State  Vector  Mode» 


Figure  C.17  Case  IV:  Estimation  Results 


x  -|  o-3  Relative  Spacecraft  Position  Error 
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Figure  C.18  Case  IV:  True  Error  of  Relative  Position  Estimates  (12-State) 


All  three  cases  satisfied  the  covariance  criterion  at  approximately  iteration 
10.  The  standard  deviation  criterion  was  not  satisfied  for  any  of  the  three  trials. 
Interestingly,  the  aresici  for  all  three  data  sets,  and  for  all  three  trials,  remained  ap¬ 
proximately  constant  to  4  significant  digits  near  iteration  15  till  program  termination 

at  imax  =  200.  The  similar  estimated  values  for  the  rotational  states,  regardless  of 
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Relative  Spacecraft  Velocity  Error 


Figure  C.19  Case  IV:  True  Error  of  Relative  Velocity  Estimates  (12-State) 


Time  [s] 

Figure  C.20  Case  IV:  True  Error  of  Angular  Velocity  Estimates  (12-State) 

the  MOI,  might  suggest  a  multiple  solution.  It  is  also  possible  that  the  differential 
equations  for  the  MOI  components  are  the  cause.  Closer  investigation  should  be  ac¬ 
complished  to  determine  the  cause  of  these  solutions  which  seem  to  be  reaching  some 
limit  of  accuracy  with  respect  to  the  residuals.  Completing  a  RSS  surface  plot  for 
this  case  might  provide  insight  to  this  unpredicted  behavior.  It  would  be  interesting 
to  use  data  sets  from  other  points  to  see  if  the  apparent  aresid  limits  changes.  It 
should  be  mentioned  that  different  values  of  k\  were  tried  with  no  significant  change 
in  the  outcome. 
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Euler  Angle  Orientation  Error 
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Figure  C.21  Case  IV:  True  Error  of  Euler  Angle  Estimates  (12-State) 


Table  C.l  Case  IV:  Estimated  9-State  Results  at  Epoch 


Estimates 

Variable 

Trial  1 

Trial  2 

Trial  3 

Truth 

[10%  MOI  Error] 

[20%  MOI  Error] 

[30%  MOI  Error] 

CJl 

rad/ s] 

0.0100 

0.0100 

0.0100 

0.01 

UJ2 

rad/ s] 

0.0001 

0.0000 

0.0000 

0 

U  3 

rad/ s] 

-0.0001 

-0.0001 

0.0000 

.001 

9i 

[rad] 

0.0004 

0.0004 

0.0000 

0 

e2 

[rad] 

-0.0107 

-0.0107 

-0.0091 

0 

0s 

[rad] 

-0.0100 

0.0007 

0.0005 

0 

A 

[ m2kg ] 

46.68 

54.55 

58.07 

45 

B 

[m2kg] 

40.62 

34.99 

40.01 

30 

C 

[m2kg] 

4.09 

12.50 

9.08 

15 

O reside 

n=i  M 

0.257 

0.257 

0.257 

0.248 

Preside 

n=2  M 

0.323 

0.327 

0.327 

0.247 

G  reside 

n=3  M 

0.287 

0.285 

0.285 

0.250 

iterations 

200 

200 

200 

- 

processing  time  [s] 

2535 

2436 

2659 

- 
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